* This is a sample main program for twpbvp, a fortran 77 package * for solving two-point boundary value problems, written by * J.R Cash and M.H. Wright. The code is fully documented * in the user's guide, available as a LaTeX file from netlib. * Machine precision is provided automatically via netlib * (in the routine d1mach). * Questions should be directed to Professor Jeff Cash, Department * of Mathematics, Imperial College, London SW7 2BZ, United Kingdom, * email j.cash@ic.ac.uk. implicit double precision (a-h,o-z) dimension fixpnt(2), ltol(2), tol(2) dimension u(2, 1000), xx(1000) dimension wrk(10000), iwrk(6000) character*30 pname logical linear, giveu, givmsh external fsub, dfsub, gsub, dgsub common/counts/nfcall * The value of eps represents the parameter epsilon * appearing in the test problem. common/probs/eps, pi logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 * double precision datan parameter (one = 1.0d+0, zero = 0.0d+0, ten = 10.0d+0) nfcall = 0 pi = 4.0d+0*datan(one) * pname = problem name pname = 'Test problem with eps = .01' eps = 1.0d-2 nudim = 2 lwrkfl = 10000 lwrkin = 6000 ncomp = 2 nlbc = 1 ntol = 1 ltol(1) = 1 tol(1) = 0.01 nfxpnt = 0 aleft = zero aright = one linear = .false. giveu = .false. givmsh = .false. write(6,771) pname write(6,772) ncomp, nlbc, ntol write(6,773) aleft, aright if (linear) write (6,774) if (.not. linear) write(6,775) if (nfxpnt .gt. 0) write(6,776) nfxpnt write(6,777) (ltol(i), i=1,ntol) write(6,778) (tol(i), i=1,ntol) call twpbvp(ncomp, nlbc, aleft, aright, * nfxpnt, fixpnt, ntol, ltol, tol, * linear, givmsh, giveu, nmsh, * xx, nudim, u, nmax, * lwrkfl, wrk, lwrkin, iwrk, * fsub, dfsub, gsub, dgsub, iflbvp) if (iflbvp .eq. 0) then write(6,901) nmsh write(6,902) if (nmsh .lt. 25) then iminc = 1 elseif (nmsh .lt. 75) then iminc = 5 elseif (nmsh .lt. 200) then iminc = 10 elseif (nmsh .lt. 1000) then iminc = 20 else iminc = 50 endif do 100 i = 1, nmsh-1, iminc write(6,900) i, xx(i), (u(j,i), j=1,ncomp) 100 continue write(6,900) nmsh, xx(nmsh), (u(j,nmsh), j=1,ncomp) endif if(iflbvp .eq. -1) then write(6,903) endif stop 771 format(1h ,a30) 772 format(1h ,'ncomp=',i5,5x,'nlbc=',i5,5x,'ntol=',i5) 773 format(1h ,'aleft =',1pe11.3,5x,'aright=',1pe11.3) 774 format(1h ,'solved as a linear problem') 775 format(1h ,'solved as a nonlinear problem') 776 format(1h ,'nfxpnt =',i5) 777 format(1h ,'ltol',5x,5i5) 778 format(1h ,'tolerances',5(1pe11.3)) 900 format(1h ,i4,5(1pe14.6)) 901 format(/, 1h ,'final solution, nmsh =',i8,/) 902 format(1h ,3x,'i',2x,'mesh point',8x,'u(1)',10x,'u(2)') 903 format(1h ,'error in the input parameters') end subroutine fsub(ncomp, x, u, f) * part of sample test function for code twpbvp, for two-point * boundary value problems. This example is documented in the * user's guide for twpbvp, a LaTeX file available from netlib. * It is called `nonlinear problem 1'. implicit double precision (a-h,o-z) dimension u(*), f(*) common/probs/eps, pi parameter (half = 0.5d+0, two = 2.0d+0) * double precision exp, sin common/counts/nfcall nfcall = nfcall + 1 * nonlinear problem 1 f(1) = u(2) f(2) = (-exp(u(1))*u(2) + * half*pi*sin(half*pi*x)*exp(two*u(1)))/eps return end subroutine dfsub(ncomp, x, u, df) * part of sample test function for code twpbvp, for two-point * boundary value problems. This example is documented in the * user's guide for twpbvp, a LaTeX file available from netlib. implicit double precision (a-h,o-z) dimension u(*), df(ncomp,*) common/probs/eps, pi parameter (half = 0.5d+0, one = 1.0d+0, zero = 0.0d+0) parameter (two = 2.0d+0) * double precision exp, sin df(1,1) = zero df(1,2) = one df(2,1) = -(exp(u(1))*u(2) * - pi*sin(pi*half*x)*exp(two*u(1)))/eps df(2,2) = -exp(u(1))/eps return end subroutine gsub(i, ncomp, u, g) * part of sample test function for code twpbvp, for two-point * boundary value problems. This example is documented in the * user's guide for twpbvp, a LaTeX file available from netlib. implicit double precision (a-h,o-z) dimension u(*) g = u(1) return end subroutine dgsub(i, ncomp, u, dg) * part of sample test function for code twpbvp, for two-point * boundary value problems. This example is documented in the * user's guide for twpbvp, a LaTeX file available from netlib. implicit double precision (a-h,o-z) dimension u(*), dg(*) parameter (one = 1.0d+0, zero = 0.0d+0) dg(1) = one dg(2) = zero return end c This is the code that has gone into netlib as a replacement. DOUBLE PRECISION FUNCTION D1MACH(I) INTEGER I C C DOUBLE-PRECISION MACHINE CONSTANTS C C D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. C C D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. C C D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. C C D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. C C D1MACH( 5) = LOG10(B) C C TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, C THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY C REMOVING THE C FROM COLUMN 1. C ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED. C (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.) C C FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), ONE OF THE FIRST C TWO SETS OF CONSTANTS BELOW SHOULD BE APPROPRIATE. IF YOU DO NOT C KNOW WHICH SET TO USE, TRY BOTH AND SEE WHICH GIVES PLAUSIBLE C VALUES. C C WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED C TO SPECIFY THE CONSTANTS EXACTLY. SOMETIMES THIS REQUIRES USING C EQUIVALENT INTEGER ARRAYS. IF YOUR COMPILER USES HALF-WORD C INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO C CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER C TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS. C C COMMENTS JUST BEFORE THE END STATEMENT (LINES STARTING WITH *) C GIVE C SOURCE FOR D1MACH. C INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) INTEGER SC, CRAY1(38), J COMMON /D9MACH/ CRAY1 C/6S C/7S SAVE SMALL, LARGE, RIGHT, DIVER, LOG10, SC C/ DOUBLE PRECISION DMACH(5) C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR BIG-ENDIAN IEEE ARITHMETIC (BINARY FORMAT) C MACHINES IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST, C SUCH AS THE AT&T 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. C SUN 3), AND MACHINES THAT USE SPARC, HP, OR IBM RISC CHIPS. C C DATA SMALL(1),SMALL(2) / 1048576, 0 / C DATA LARGE(1),LARGE(2) / 2146435071, -1 / C DATA RIGHT(1),RIGHT(2) / 1017118720, 0 / C DATA DIVER(1),DIVER(2) / 1018167296, 0 / C DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 /, SC/987/ C C MACHINE CONSTANTS FOR LITTLE-ENDIAN (BINARY) IEEE ARITHMETIC C MACHINES IN WHICH THE LEAST SIGNIFICANT BYTE IS STORED FIRST, C E.G. IBM PCS AND OTHER MACHINES THAT USE INTEL 80X87 OR DEC C ALPHA CHIPS. C C DATA SMALL(1),SMALL(2) / 0, 1048576 / C DATA LARGE(1),LARGE(2) / -1, 2146435071 / C DATA RIGHT(1),RIGHT(2) / 0, 1017118720 / C DATA DIVER(1),DIVER(2) / 0, 1018167296 / C DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 /, SC/987/ C C MACHINE CONSTANTS FOR AMDAHL MACHINES. C C DATA SMALL(1),SMALL(2) / 1048576, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 856686592, 0 / C DATA DIVER(1),DIVER(2) / 873463808, 0 / C DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 /, SC/987/ C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 /, SC/987/ C C MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / 00564000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C C DATA LARGE(1) / 37757777777777777777B / C DATA LARGE(2) / 37157777777777777774B / C C DATA RIGHT(1) / 15624000000000000000B / C DATA RIGHT(2) / 00000000000000000000B / C C DATA DIVER(1) / 15634000000000000000B / C DATA DIVER(2) / 00000000000000000000B / C C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B /, SC/987/ C C MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. C C DATA SMALL(1) / O"00564000000000000000" / C DATA SMALL(2) / O"00000000000000000000" / C C DATA LARGE(1) / O"37757777777777777777" / C DATA LARGE(2) / O"37157777777777777774" / C C DATA RIGHT(1) / O"15624000000000000000" / C DATA RIGHT(2) / O"00000000000000000000" / C C DATA DIVER(1) / O"15634000000000000000" / C DATA DIVER(2) / O"00000000000000000000" / C C DATA LOG10(1) / O"17164642023241175717" / C DATA LOG10(2) / O"16367571421742254654" /, SC/987/ C C MACHINE CONSTANTS FOR CONVEX C-1 C C DATA SMALL(1),SMALL(2) / '00100000'X, '00000000'X / C DATA LARGE(1),LARGE(2) / '7FFFFFFF'X, 'FFFFFFFF'X / C DATA RIGHT(1),RIGHT(2) / '3CC00000'X, '00000000'X / C DATA DIVER(1),DIVER(2) / '3CD00000'X, '00000000'X / C DATA LOG10(1),LOG10(2) / '3FF34413'X, '509F79FF'X /, SC/987/ C C MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. C C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777776B / C C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B /, SC/987/ C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C C SMALL, LARGE, RIGHT, DIVER, LOG10 SHOULD BE DECLARED C INTEGER SMALL(4), LARGE(4), RIGHT(4), DIVER(4), LOG10(4) C C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - C STATIC DMACH(5) C C DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ C DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ C DATA LOG10/40423K,42023K,50237K,74776K/, SC/987/ C C MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 C C DATA SMALL(1),SMALL(2) / '20000000, '00000201 / C DATA LARGE(1),LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1),RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1),DIVER(2) / '20000000, '00000334 / C DATA LOG10(1),LOG10(2) / '23210115, '10237777 /, SC/987/ C C MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. C C DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /, SC/987/ C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. C C DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF /, SC/987/ C C MACHINE CONSTANTS FOR THE INTERDATA 8/32 C WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. C C FOR THE INTERDATA FORTRAN VII COMPILER REPLACE C THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. C C DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000' / C DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000' / C DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF' /, SC/987/ C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). C C DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 /, SC/987/ C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). C C DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1),LOG10(2) / "177464202324, "047674776746 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1),SMALL(2) / 8388608, 0 / C DATA LARGE(1),LARGE(2) / 2147483647, -1 / C DATA RIGHT(1),RIGHT(2) / 612368384, 0 / C DATA DIVER(1),DIVER(2) / 620756992, 0 / C DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /, SC/987/ C C DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 /, SC/987/ C C MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C SMALL, LARGE, RIGHT, DIVER, LOG10 SHOULD BE DECLARED C INTEGER SMALL(4), LARGE(4), RIGHT(4), DIVER(4), LOG10(4) C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA SMALL(3),SMALL(4) / 0, 0 / C C DATA LARGE(1),LARGE(2) / 32767, -1 / C DATA LARGE(3),LARGE(4) / -1, -1 / C C DATA RIGHT(1),RIGHT(2) / 9344, 0 / C DATA RIGHT(3),RIGHT(4) / 0, 0 / C C DATA DIVER(1),DIVER(2) / 9472, 0 / C DATA DIVER(3),DIVER(4) / 0, 0 / C C DATA LOG10(1),LOG10(2) / 16282, 8346 / C DATA LOG10(3),LOG10(4) / -31493, -12296 /, SC/987/ C C DATA SMALL(1),SMALL(2) / O000200, O000000 / C DATA SMALL(3),SMALL(4) / O000000, O000000 / C C DATA LARGE(1),LARGE(2) / O077777, O177777 / C DATA LARGE(3),LARGE(4) / O177777, O177777 / C C DATA RIGHT(1),RIGHT(2) / O022200, O000000 / C DATA RIGHT(3),RIGHT(4) / O000000, O000000 / C C DATA DIVER(1),DIVER(2) / O022400, O000000 / C DATA DIVER(3),DIVER(4) / O000000, O000000 / C C DATA LOG10(1),LOG10(2) / O037632, O020232 / C DATA LOG10(3),LOG10(4) / O102373, O147770 /, SC/987/ C C MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS C WITH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS, C SUPPLIED BY IGOR BRAY. C C DATA SMALL(1),SMALL(2) / :10000000000, :00000100001 / C DATA LARGE(1),LARGE(2) / :17777777777, :37777677775 / C DATA RIGHT(1),RIGHT(2) / :10000000000, :00000000122 / C DATA DIVER(1),DIVER(2) / :10000000000, :00000000123 / C DATA LOG10(1),LOG10(2) / :11504046501, :07674600177 /, SC/987/ C C MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000 C C DATA SMALL(1),SMALL(2) / $00000000, $00100000 / C DATA LARGE(1),LARGE(2) / $FFFFFFFF, $7FEFFFFF / C DATA RIGHT(1),RIGHT(2) / $00000000, $3CA00000 / C DATA DIVER(1),DIVER(2) / $00000000, $3CB00000 / C DATA LOG10(1),LOG10(2) / $509F79FF, $3FD34413 /, SC/987/ C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. C C DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX UNIX F77 COMPILER C C DATA SMALL(1),SMALL(2) / 128, 0 / C DATA LARGE(1),LARGE(2) / -32769, -1 / C DATA RIGHT(1),RIGHT(2) / 9344, 0 / C DATA DIVER(1),DIVER(2) / 9472, 0 / C DATA LOG10(1),LOG10(2) / 546979738, -805796613 /, SC/987/ C C MACHINE CONSTANTS FOR THE VAX-11 WITH C FORTRAN IV-PLUS COMPILER C C DATA SMALL(1),SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1),LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1),RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1),DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1),LOG10(2) / Z209A3F9A, ZCFF884FB /, SC/987/ C C MACHINE CONSTANTS FOR VAX/VMS VERSION 2.2 C C DATA SMALL(1),SMALL(2) / '80'X, '0'X / C DATA LARGE(1),LARGE(2) / 'FFFF7FFF'X, 'FFFFFFFF'X / C DATA RIGHT(1),RIGHT(2) / '2480'X, '0'X / C DATA DIVER(1),DIVER(2) / '2500'X, '0'X / C DATA LOG10(1),LOG10(2) / '209A3F9A'X, 'CFF884FB'X /, SC/987/ C C *** ISSUE STOP 779 IF ALL DATA STATEMENTS ARE COMMENTED... IF (SC .NE. 987) THEN DMACH(1) = 1.D13 IF ( SMALL(1) .EQ. 1117925532 * .AND. SMALL(2) .EQ. -448790528) THEN * *** IEEE BIG ENDIAN *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2146435071 LARGE(2) = -1 RIGHT(1) = 1017118720 RIGHT(2) = 0 DIVER(1) = 1018167296 DIVER(2) = 0 LOG10(1) = 1070810131 LOG10(2) = 1352628735 ELSE IF ( SMALL(2) .EQ. 1117925532 * .AND. SMALL(1) .EQ. -448790528) THEN * *** IEEE LITTLE ENDIAN *** SMALL(2) = 1048576 SMALL(1) = 0 LARGE(2) = 2146435071 LARGE(1) = -1 RIGHT(2) = 1017118720 RIGHT(1) = 0 DIVER(2) = 1018167296 DIVER(1) = 0 LOG10(2) = 1070810131 LOG10(1) = 1352628735 ELSE IF ( SMALL(1) .EQ. -2065213935 * .AND. SMALL(2) .EQ. 10752) THEN * *** VAX WITH D_FLOATING *** SMALL(1) = 128 SMALL(2) = 0 LARGE(1) = -32769 LARGE(2) = -1 RIGHT(1) = 9344 RIGHT(2) = 0 DIVER(1) = 9472 DIVER(2) = 0 LOG10(1) = 546979738 LOG10(2) = -805796613 ELSE IF ( SMALL(1) .EQ. 1267827943 * .AND. SMALL(2) .EQ. 704643072) THEN * *** IBM MAINFRAME *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2147483647 LARGE(2) = -1 RIGHT(1) = 856686592 RIGHT(2) = 0 DIVER(1) = 873463808 DIVER(2) = 0 LOG10(1) = 1091781651 LOG10(2) = 1352628735 ELSE IF ( SMALL(1) .EQ. 1120022684 * .AND. SMALL(2) .EQ. -448790528) THEN * *** CONVEX C-1 *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2147483647 LARGE(2) = -1 RIGHT(1) = 1019215872 RIGHT(2) = 0 DIVER(1) = 1020264448 DIVER(2) = 0 LOG10(1) = 1072907283 LOG10(2) = 1352628735 ELSE IF ( SMALL(1) .EQ. 815547074 * .AND. SMALL(2) .EQ. 58688) THEN * *** VAX G-FLOATING *** SMALL(1) = 16 SMALL(2) = 0 LARGE(1) = -32769 LARGE(2) = -1 RIGHT(1) = 15552 RIGHT(2) = 0 DIVER(1) = 15568 DIVER(2) = 0 LOG10(1) = 1142112243 LOG10(2) = 2046775455 ELSE DMACH(2) = 1.D27 + 1 DMACH(3) = 1.D27 LARGE(2) = LARGE(2) - RIGHT(2) IF (LARGE(2) .EQ. 64 .AND. SMALL(2) .EQ. 0) THEN CRAY1(1) = 67291416 DO 10 J = 1, 20 10 CRAY1(J+1) = CRAY1(J) + CRAY1(J) CRAY1(22) = CRAY1(21) + 321322 DO 20 J = 22, 37 20 CRAY1(J+1) = CRAY1(J) + CRAY1(J) IF (CRAY1(38) .EQ. SMALL(1)) THEN * *** CRAY *** * SMALL(1) = 2332160919536140288 SMALL(1) = 2332160 SMALL(1) = 1000000*SMALL(1) + 919536 SMALL(1) = 1000000*SMALL(1) + 140288 SMALL(2) = 0 * LARGE(1) = 6917247552664371199 LARGE(1) = 6917247 LARGE(1) = 1000000*LARGE(1) + 552664 LARGE(1) = 1000000*LARGE(1) + 371199 * LARGE(2) = 281474976710654 LARGE(2) = 28147497 LARGE(2) = 10000000*LARGE(2) + 6710654 * RIGHT(1) = 4585649583081652224 RIGHT(1) = 4585649 RIGHT(1) = 1000000*RIGHT(1) + 583081 RIGHT(1) = 1000000*RIGHT(1) + 652224 RIGHT(2) = 0 * DIVER(1) = 4585931058058362880 DIVER(1) = 4585931 DIVER(1) = 1000000*DIVER(1) + 058058 DIVER(1) = 1000000*DIVER(1) + 362880 DIVER(2) = 0 * LOG10(1) = 4611574008272714703 LOG10(1) = 4611574 LOG10(1) = 1000000*LOG10(1) + 8272 LOG10(1) = 1000000*LOG10(1) + 714703 * LOG10(2) = 272234615232940 LOG10(2) = 27223461 LOG10(2) = 10000000*LOG10(2) + 5232940 ELSE WRITE(*,9000) STOP 779 END IF ELSE WRITE(*,9000) STOP 779 END IF END IF SC = 987 END IF C C *** ISSUE STOP 778 IF ALL DATA STATEMENTS ARE OBVIOUSLY WRONG... IF (DMACH(4) .GE. 1.0D0) STOP 778 *C/6S *C IF (I .LT. 1 .OR. I .GT. 5) *C 1 CALL SETERR(24HD1MACH - I OUT OF BOUNDS,24,1,2) *C/7S * IF (I .LT. 1 .OR. I .GT. 5) * 1 CALL SETERR('D1MACH - I OUT OF BOUNDS',24,1,2) *C/ IF (I .LT. 1 .OR. I .GT. 5) THEN WRITE(*,*) 'D1MACH(I): I =',I,' is out of bounds.' STOP END IF D1MACH = DMACH(I) RETURN C/6S C9000 FORMAT(/46H Adjust D1MACH by uncommenting data statements/ C *30H appropriate for your machine.) C/7S 9000 FORMAT(/' Adjust D1MACH by uncommenting data statements'/ *' appropriate for your machine.') C/ C * /* C source for D1MACH -- remove the * in column 1 */ *#include *#include *#include * *double d1mach_(long *i) *{ * switch(*i){ * case 1: return DBL_MIN; * case 2: return DBL_MAX; * case 3: return DBL_EPSILON/FLT_RADIX; * case 4: return DBL_EPSILON; * case 5: return log10(FLT_RADIX); * } * * fprintf(stderr, "invalid argument: d1mach(%ld)\n", *i); * exit(1); * return 0; /* for compilers that complain of missing return values */ * } END subroutine sprt(n, array) * sprt prints an array. implicit double precision (a-h,o-z) dimension array(n) write(6,900) (array(i), i=1,n) 900 format(1h ,(7(1pe11.3))) return end subroutine mprt(nrowd, nrow, ncol, array) * mprt prints a matrix. implicit double precision (a-h,o-z) dimension array(nrowd, ncol) do 400 i = 1, nrow write(6,900) i,(array(i,j), j=1,ncol) 400 continue 900 format(1h ,i5,(6(1pe11.3))) return end subroutine twpbvp(ncomp, nlbc, aleft, aright, * nfxpnt, fixpnt, ntol, ltol, tol, * linear, givmsh, giveu, nmsh, * xx, nudim, u, nmax, * lwrkfl, wrk, lwrkin, iwrk, * fsub, dfsub, gsub, dgsub, iflbvp) * subroutine twpbvp is intended to solve two-point boundary * value problems. For information about its use, parameters, * etc., see the user's guide, a LaTeX file available from netlib. implicit double precision (a-h,o-z) dimension fixpnt(*), ltol(*), tol(*) dimension xx(*), u(nudim,*) dimension wrk(lwrkfl), iwrk(lwrkin) logical linear, givmsh, giveu external fsub, dfsub, gsub, dgsub logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 intrinsic abs, min parameter ( zero = 0.0d+0 ) * Check for invalid input parameters. If any parameters are * invalid, exit with the flag iflbvp set to -1. iflbvp = -1 if (ncomp .le. 0) return if (nlbc .lt. 0 .or. nlbc .gt. ncomp) return if (aleft .ge. aright) return if (nfxpnt .lt. 0) return if (givmsh .and. nmsh .lt. nfxpnt+2) return if (givmsh .and. xx(1) .ne. aleft) return if (givmsh .and. xx(nmsh) .ne. aright) return if (nfxpnt .gt. 0) then if (fixpnt(1) .le. aleft) return if (fixpnt(nfxpnt) .ge. aright) return do 50 i = 1, nfxpnt-1 if (fixpnt(i+1) .le. fixpnt(i)) return 50 continue endif if (ntol .lt. 1) return do 60 i = 1, ntol if (ltol(i) .lt. 0 .or. ltol(i) .gt. ncomp) return if (tol(i) .le. zero) return 60 continue if (giveu .and. .not. givmsh) return if (nudim .le. 0) return if (lwrkfl .le. 0 .or. lwrkin .le. 0) return * Calculate maximum number of mesh points possible with the * given floating-point and integer workspace. isp = lwrkfl - 2*ntol - 9*ncomp - 4*ncomp*ncomp iden = 4*ncomp*ncomp + 12*ncomp + 3 nmax1 = isp/iden isp = lwrkin - ncomp nmax2 = isp/(ncomp+2) nmax = min(nmax1, nmax2) if (iprint .ge. 0) write(6,901) nmax 901 format(1h ,'nmax from workspace =',i8) if (nmax .le. 1) return * Partition floating point workspace. irhs = 1 lrhs = ncomp*nmax itpblk = irhs + lrhs ltpblk = ncomp*nlbc ibtblk = itpblk + ltpblk lbtblk = ncomp*(ncomp - Nlbc) iajac = ibtblk + lbtblk lajac = 2*ncomp*ncomp*nmax ibhold = iajac + lajac lbhold = ncomp*ncomp*nmax ichold = ibhold + lbhold lchold = ncomp*ncomp*nmax ifval = ichold + lchold lfval = ncomp*nmax idef = ifval + lfval ldef = ncomp*(nmax-1) idefex = idef + ldef ldefex = ncomp*(nmax-1) * def6 uses the same space as defexp idef6 = idefex ldef6 = ncomp*(nmax-1) idefim = idef6 + ldef6 ldefim = ncomp*(nmax-1) * def8 uses the same space as defimp idef8 = idefim ldef8 = ncomp*(nmax-1) iusve = idef8 + ldef8 lusve = ncomp*nmax iuold = iusve + lusve luold = ncomp*nmax itmrhs = iuold + luold ltmrhs = ncomp*nmax irhtri = itmrhs + ltmrhs lrhtri = ncomp*nmax idelu = irhtri + lrhtri ldelu = ncomp*nmax ixmer = idelu + ldelu lxmer = ncomp*nmax * rerr occupies the same space as xmerit irerr = ixmer lrerr = ncomp*nmax iutri = irerr + lrerr lutri = ncomp*nmax iermx = iutri + lutri lermx = nmax irtdc = iermx + lermx lrtdc = nmax ixxold = irtdc + lrtdc lxxold = nmax iuint = ixxold + lxxold luint = ncomp iftmp = iuint + luint lftmp = ncomp idgtm = iftmp + lftmp ldgtm = ncomp idftm1 = idgtm + ldgtm ldftm1 = ncomp*ncomp idftm2 = idftm1 + ldftm1 ldftm2 = ncomp*ncomp itmp = idftm2 + ldftm2 ltmp = ncomp*8 idsq = itmp + ltmp ldsq = ncomp*ncomp idexr = idsq + ldsq ldexr = ncomp ietst6 = idexr + ldexr letst6 = ntol ietst8 = ietst6 + letst6 letst8 = ntol ilast = ietst8 + letst8 if (iprint .eq. 1) write(6,903) ilast 903 format(1h ,'ilast',i10) * Partition integer workspace. iiref = 1 liref = nmax iihcom = iiref + liref lihcom = nmax iipvbk = iihcom + lihcom lipvbk = ncomp*nmax iipvlu = iipvbk + lipvbk lipvlu = ncomp call bvpsol(ncomp, nmsh, nlbc, aleft, aright, * nfxpnt, fixpnt, ntol, ltol, tol, nmax, linear, * giveu, givmsh, xx, nudim, u, * wrk(idefex), wrk(idefim), wrk(idef), wrk(idelu), * wrk(irhs), wrk(ifval), * wrk(itpblk), wrk(ibtblk), wrk(iajac), wrk(ibhold), * wrk(ichold), iwrk(iipvbk), iwrk(iipvlu), * wrk(iuint), wrk(iftmp), wrk(itmrhs), * wrk(idftm1), wrk(idftm2), wrk(idgtm), * wrk(iutri), wrk(irhtri), wrk(ixmer), * wrk(ixxold), wrk(iuold), wrk(iusve), * wrk(itmp), wrk(idsq), wrk(idexr), wrk(irtdc), * wrk(irerr), wrk(ietst6), wrk(ietst8), wrk(iermx), * iwrk(iihcom), iwrk(iiref), wrk(idef6), wrk(idef8), * fsub,dfsub,gsub,dgsub,iflbvp) return end C----------------------------------------------------------------------- C C C----------------------------------------------------------------------- C SUBROUTINE COLROW (N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK, * NBLOKS,BOTBLK,NRWBOT,PIVOT,B,X,IFLAG) C C*************************************************************** C C THIS PROGRAM SOLVES THE LINEAR SYSTEM A*X = B WHERE A IS C AN ALMOST BLOCK DIAGONAL MATRIX OF THE FORM C C TOPBLK C ARRAY(1) C ARRAY(2) C . C . C . C . C ARRAY(NBLOKS) C BOTBLK C C WHERE C TOPBLK IS NRWTOP BY NOVRLP C ARRAY(K), K=1,NBLOKS, ARE NRWBLK BY NRWBLK+NOVRLP C BOTBLK IS NRWBOT BY NOVRLP, C AND C NOVRLP = NRWTOP + NRWBOT C WITH C NOVRLP.LE.NRWBLK . C C THE LINEAR SYSTEM IS OF ORDER N = NBLOKS*NRWBLK + NOVRLP. C C THE METHOD IMPLEMENTED IS BASED ON GAUSS ELIMINATION WITH C ALTERNATE ROW AND COLUMN ELIMINATION WITH PARTIAL PIVOTING, C WHICH PRODUCES A STABLE DECOMPOSITION OF THE MATRIX A C WITHOUT INTRODUCING FILL-IN. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C TO OBTAIN A SINGLE PRECISION VERSION OF THIS PACKAGE, REMOVE C ALL DOUBLE PRECISION STATEMENTS. THERE IS ONE SUCH STATEMENT C IN C O L R O W, THREE IN C R D C M P, AND TWO IN C R S O L V. C IN ADDITION, REFERENCES TO BUILT-IN FUNCTIONS DABS AND DMAX1 C MUST BE REPLACED BY DABS AND DMAX1, RESPECTIVELY. ABS OCCURS C NINE TIMES, IN C R D C M P. DMAX1 OCCURS FOUR TIMES, IN C C R D C M P. FINALLY, ZERO IS INITIALISED TO 0.D+0 IN A C DATA STATEMENT IN C R D C M P. THIS MUST BE REPLACED BY: C DATA ZERO/0.0/ C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** PARAMETERS ***** C C *** ON ENTRY ... C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM, C GIVEN BY NBLOKS*NRWBLK + NOVRLP C C TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP) C THE FIRST BLOCK OF THE ALMOST BLOCK C DIAGONAL MATRIX A C C NRWTOP - INTEGER C NUMBER OF ROWS IN THE BLOCK TOPBLK C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN WHICH SUCC- C ESSIVE BLOCKS OVERLAP, WHERE C NOVRLP = NRWTOP + NRWBOT C C ARRAY - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS) C ARRAY(,K) CONTAINS THE K-TH NRWBLK C BY NCLBLK BLOCK OF THE MATRIX A C C NRWBLK - INTEGER C NUMBER OF ROWS IN K-TH BLOCK C C NCLBLK - INTEGER C NUMBER OF COLUMNS IN K-TH BLOCK C C NBLOKS - INTEGER C NUMBER OF NRWBLK BY NCLBLK BLOCKS IN C THE MATRIX A C C BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP) C THE LAST BLOCK OF THE MATRIX A C C NRWBOT - INTEGER C NUMBER OF ROWS IN THE BLOCK BOTBLK C C PIVOT - INTEGER(N) C WORK SPACE C C B - DOUBLE PRECISION(N) C THE RIGHT HAND SIDE VECTOR C C X - DOUBLE PRECISION(N) C WORK SPACE C C *** ON RETURN ... C C TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE C DESIRED DECOMPOSITION OF THE MATRIX A C (IF IFLAG = 0) C C PIVOT - INTEGER(N) C RECORDS THE PIVOTING INDICES DETER- C MINED IN THE DECOMPOSITION C C X - DOUBLE PRECISION(N) C THE SOLUTION VECTOR (IF IFLAG = 0) C C IFLAG - INTEGER C = 1, IF INPUT PARAMETERS ARE INVALID C = -1, IF MATRIX IS SINGULAR C = 0, OTHERWISE C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** AUXILIARY PROGRAMS ***** C C CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS, C * BOTBLK,NRWBOT,PIVOT,IFLAG) C - DECOMPOSES THE MATRIX A USING MODIFIED C ALTERNATE ROW AND COLUMN ELIMINATON WITH C PARTIAL PIVOTING, AND IS USED FOR THIS C PURPOSE IN C O L R O W. C THE ARGUMENTS ARE AS IN C O L R O W. C C CRSLVE(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS, C * BOTBLK,NRWBOT,PIVOT,B,X) C - SOLVES THE SYSTEM A*X = B ONCE A IS DECOMPOSED. C THE ARGUMENTS ARE ALLAS IN C O L R O W. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C THE SUBROUTINE C O L R O W AUTOMATICALLY SOLVES THE C INPUT SYSTEM WHEN IFLAG=0. C O L R O W IS CALLED ONLY ONCE C FOR A GIVEN SYSTEM. THE SOLUTION FOR A SEQUENCE OF P RIGHT C HAND SIDES CAN BE OBTAINED BY ONE CALL TO C O L R O W AND C P-1 CALLS TO CRSLVE ONLY. SINCE THE ARRAYS TOPBLK,ARRAY, C BOTBLK AND PIVOT CONTAIN THE DECOMPOSITION OF THE GIVEN C COEFFICIENT MATRIX AND PIVOTING INFORMATION ON RETURN FROM C C O L R O W , THEY MUST NOT BE ALTERED BETWEEN SUCCESSIVE C CALLS TO CRSLVE WITH THE SAME LEFT HAND SIDES. FOR THE C SAME REASON, IF THE USER WISHES TO SAVE THE COEFFICIENT C MATRIX, THE ARRAYS TOPBLK,ARRAY,BOTBLK MUST BE COPIED C BEFORE A CALL TO C O L R O W . C C*********************************************************************** C*************************************************************** C IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER PIVOT(1) DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1),BOTBLK(NRWBOT,1) * ,B(1),X(1) CALL CRDCMP (N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS, * BOTBLK,NRWBOT,PIVOT,IFLAG) IF (IFLAG.NE.0) RETURN CALL CRSLVE (TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS, * BOTBLK,NRWBOT,PIVOT,B,X) RETURN END SUBROUTINE CRDCMP(N,TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK, * NBLOKS,BOTBLK,NRWBOT,PIVOT,IFLAG) C C*************************************************************** C C C R D C M P DECOMPOSES THE ALMOST BLOCK DIAGONAL MATRIX A C USING MODIFIED ALTERNATE ROW AND COLUMN ELIMINATION WITH C PARTIAL PIVOTING. THE MATRIX A IS STORED IN THE ARRAYS C TOPBLK, ARRAY, AND BOTBLK. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** PARAMETERS ***** C C *** ON ENTRY ... C C N - INTEGER C THE ORDER OF THE LINEAR SYSTEM, C GIVEN BY NBLOKS*NRWBLK + NOVRLP C C TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP) C THE FIRST BLOCK OF THE ALMOST BLOCK C DIAGONAL MATRIX A TO BE DECOMPOSED C C NRWTOP - INTEGER C NUMBER OF ROWS IN THE BLOCK TOPBLK C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN WHICH SUCC- C ESSIVE BLOCKS OVERLAP, WHERE C NOVRLP = NRWTOP + NRWBOT C C ARRAY - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS) C ARRAY(,,K) CONTAINS THE K-TH NRWBLK C BY NCLBLK BLOCK OF THE MATRIX A C C NRWBLK - INTEGER C NUMBER OF ROWS IN K-TH BLOCK C C NCLBLK - INTEGER C NUMBER OF COLUMNS IN K-TH BLOCK C C NBLOKS - INTEGER C NUMBER OF NRWBLK BY NCLBLK BLOCKS IN C THE MATRIX A C C BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP) C THE LAST BLOCK OF THE MATRIX A C C NRWBOT - INTEGER C NUMBER OF ROWS IN THE BLOCK BOTBLK C C PIVOT - INTEGER(N) C WORK SPACE C C *** ON RETURN ... C C TOPBLK,ARRAY,BOTBLK - ARRAYS CONTAINING THE C DESIRED DECOMPOSITION OF THE MATRIX A C (IF IFLAG = 0) C C PIVOT - INTEGER(N) C RECORDS THE PIVOTING INDICES DETER- C MINED IN THE DECOMPOSITION C C IFLAG - INTEGER C = 1, IF INPUT PARAMETERS ARE INVALID C = -1, IF MATRIX IS SINGULAR C = 0, OTHERWISE C C*************************************************************** C IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER PIVOT(1) DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1),BOTBLK(NRWBOT,1) DATA ZERO / 0.0D+0 / C C*************************************************************** C C **** DEFINE THE CONSDTANTS USED THROUGHOUT **** C C*************************************************************** C IFLAG = 0 PIVMAX = ZERO NRWTP1 = NRWTOP+1 NROWEL = NRWBLK-NRWTOP NRWEL1 = NROWEL+1 NVRLP0 = NOVRLP-1 C C*************************************************************** C C **** CHECK VALIDITY OF THE INPUT PARAMETERS.... C C IF PARAMETERS ARE INVALID THEN TERMINATE AT 10; C ELSE CONTINUE AT 100. C C*************************************************************** C IF (N.NE.NBLOKS*NRWBLK+NOVRLP) GO TO 10 IF (NOVRLP.NE.NRWTOP+NRWBOT) GO TO 10 IF (NCLBLK.NE.NOVRLP+NRWBLK) GO TO 10 IF (NOVRLP.GT.NRWBLK) GO TO 10 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C PARAMETERS ARE ACCEPTABLE - CONTINUE AT 100. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C GO TO 20 10 CONTINUE C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C PARAMETERS ARE INVALID. SET IFLAG = 1, AND TERMINATE C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IFLAG = 1 RETURN 20 CONTINUE C C*************************************************************** C C **** FIRST, IN TOPBLK.... C C*************************************************************** C C *** APPLY NRWTOP COLUMN ELIMINATIONS WITH COLUMN C PIVOTING .... C C*************************************************************** C DO 110 I = 1, NRWTOP IPLUS1 = I+1 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DETERMINE COLUMN PIVOT AND PIVOT INDEX C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IPVT = I COLMAX = DABS(TOPBLK(I,I)) DO 30 J = IPLUS1, NOVRLP TEMPIV = DABS(TOPBLK(I,J)) IF (TEMPIV.LE.COLMAX) GO TO 30 IPVT = J COLMAX = TEMPIV 30 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C TEST FOR SINGULARITY: C C IF SINGULAR THEN TERMINATE AT 1000; C ELSE CONTINUE. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF (PIVMAX+COLMAX.EQ.PIVMAX) GO TO 410 PIVMAX = DMAX1(COLMAX,PIVMAX) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IF NECESSARY INTERCHANGE COLUMNS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PIVOT(I) = IPVT IF (IPVT.EQ.I) GO TO 60 DO 40 L = I, NRWTOP SWAP = TOPBLK(L,IPVT) TOPBLK(L,IPVT) = TOPBLK(L,I) TOPBLK(L,I) = SWAP 40 CONTINUE DO 50 L = 1, NRWBLK SWAP = ARRAY(L,IPVT,1) ARRAY(L,IPVT,1) = ARRAY(L,I,1) ARRAY(L,I,1) = SWAP 50 CONTINUE 60 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C COMPUTE MULTIPLIERS AND PERFORM COLUMN C ELIMINATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COLPIV = TOPBLK(I,I) DO 100 J = IPLUS1, NOVRLP COLMLT = TOPBLK(I,J)/COLPIV TOPBLK(I,J) = COLMLT IF (IPLUS1.GT.NRWTOP) GO TO 80 DO 70 L = IPLUS1, NRWTOP TOPBLK(L,J) = TOPBLK(L,J)-COLMLT*TOPBLK(L,I) 70 CONTINUE 80 CONTINUE DO 90 L = 1, NRWBLK ARRAY(L,J,1) = ARRAY(L,J,1)-COLMLT*ARRAY(L,I,1) 90 CONTINUE 100 CONTINUE 110 CONTINUE C C*************************************************************** C C **** IN EACH BLOCK ARRAY(,,K).... C C*************************************************************** C INCR = 0 DO 320 K = 1, NBLOKS KPLUS1 = K+1 C C ***************************************************** C C *** FIRST APPLY NRWBLK-NRWTOP ROW ELIMINATIONS WITH C ROW PIVOTING.... C C ***************************************************** C DO 180 J = NRWTP1, NRWBLK JPLUS1 = J+1 JMINN = J-NRWTOP C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DETERMINE ROW PIVOT AND PIVOT INDEX C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IPVT = JMINN ROWMAX = DABS(ARRAY(JMINN,J,K)) LOOP = JMINN+1 DO 120 I = LOOP, NRWBLK TEMPIV = DABS(ARRAY(I,J,K)) IF (TEMPIV.LE.ROWMAX) GO TO 120 IPVT = I ROWMAX = TEMPIV 120 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C TEST FOR SINGULARITY: C C IF SINGULAR THEN TERMINATE AT 1000; C ELSE CONTINUE. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF (PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 410 PIVMAX = DMAX1(ROWMAX,PIVMAX) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IF NECESSARY INTERCHANGE ROWS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INCRJ = INCR+J PIVOT(INCRJ) = INCR+IPVT+NRWTOP IF (IPVT.EQ.JMINN) GO TO 140 DO 130 L = J, NCLBLK SWAP = ARRAY(IPVT,L,K) ARRAY(IPVT,L,K) = ARRAY(JMINN,L,K) ARRAY(JMINN,L,K) = SWAP 130 CONTINUE 140 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C COMPUTE MULTIPLERS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ROWPIV = ARRAY(JMINN,J,K) DO 150 I = LOOP, NRWBLK ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV 150 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C PERFORM ROW ELIMINATION WITH COLUMN INDEXING C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 170 L = JPLUS1, NCLBLK ROWMLT = ARRAY(JMINN,L,K) DO 160 I = LOOP, NRWBLK ARRAY(I,L,K) = ARRAY(I,L,K)-ROWMLT*ARRAY(I,J,K) 160 CONTINUE 170 CONTINUE 180 CONTINUE C C ***************************************************** C C *** NOW APPLY NRWTOP COLUMN ELIMINATIONS WITH C COLUMN PIVOTING.... C C ***************************************************** C DO 310 I = NRWEL1, NRWBLK IPLUSN = I+NRWTOP IPLUS1 = I+1 C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DETERMINE COLUMN PIVOT AND PIVOT INDEX C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IPVT = IPLUSN COLMAX = DABS(ARRAY(I,IPVT,K)) LOOP = IPLUSN+1 DO 190 J = LOOP, NCLBLK TEMPIV = DABS(ARRAY(I,J,K)) IF (TEMPIV.LE.COLMAX) GO TO 190 IPVT = J COLMAX = TEMPIV 190 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C TEST FOR SINGULARITY: C C IF SINGULAR THEN TERMINATE AT 1000; C ELSE CONTINUE. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF (PIVMAX+COLMAX.EQ.PIVMAX) GO TO 410 PIVMAX = DMAX1(COLMAX,PIVMAX) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IF NECESSARY INTERCHANGE COLUMNS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INCRN = INCR+IPLUSN PIVOT(INCRN) = INCR+IPVT IRWBLK = IPLUSN-NRWBLK IF (IPVT.EQ.IPLUSN) GO TO 240 DO 200 L = I, NRWBLK SWAP = ARRAY(L,IPVT,K) ARRAY(L,IPVT,K) = ARRAY(L,IPLUSN,K) ARRAY(L,IPLUSN,K) = SWAP 200 CONTINUE IPVBLK = IPVT-NRWBLK IF (K.EQ.NBLOKS) GO TO 220 DO 210 L = 1, NRWBLK SWAP = ARRAY(L,IPVBLK,KPLUS1) ARRAY(L,IPVBLK,KPLUS1) = ARRAY(L,IRWBLK,KPLUS1) ARRAY(L,IRWBLK,KPLUS1) = SWAP 210 CONTINUE GO TO 240 220 CONTINUE DO 230 L = 1, NRWBOT SWAP = BOTBLK(L,IPVBLK) BOTBLK(L,IPVBLK) = BOTBLK(L,IRWBLK) BOTBLK(L,IRWBLK) = SWAP 230 CONTINUE 240 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C COMPUTE MULTIPLIERS AND PERFORM COLUMN C ELIMINATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C COLPIV = ARRAY(I,IPLUSN,K) DO 300 J = LOOP, NCLBLK COLMLT = ARRAY(I,J,K)/COLPIV ARRAY(I,J,K) = COLMLT IF (I.EQ.NRWBLK) GO TO 260 DO 250 L = IPLUS1, NRWBLK ARRAY(L,J,K) = ARRAY(L,J,K)-COLMLT*ARRAY(L,IPLUSN,K) 250 CONTINUE 260 CONTINUE JRWBLK = J-NRWBLK IF (K.EQ.NBLOKS) GO TO 280 DO 270 L = 1, NRWBLK ARRAY(L,JRWBLK,KPLUS1) = ARRAY(L,JRWBLK,KPLUS1)-COLMLT * *ARRAY(L,IRWBLK,KPLUS1) 270 CONTINUE GO TO 300 280 CONTINUE DO 290 L = 1, NRWBOT BOTBLK(L,JRWBLK) = BOTBLK(L,JRWBLK)-COLMLT*BOTBLK(L, * IRWBLK) 290 CONTINUE 300 CONTINUE 310 CONTINUE INCR = INCR+NRWBLK 320 CONTINUE C C*************************************************************** C C **** FINALLY, IN BOTBLK.... C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C *** APPLY NRWBOT ROW ELIMINATIONS WITH ROW C PIVOTING.... C C IF BOT HAS JUST ONE ROW GO TO 500 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF (NRWBOT.EQ.1) GO TO 400 DO 390 J = NRWTP1, NVRLP0 JPLUS1 = J+1 JMINN = J-NRWTOP C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DETERMINE ROW PIVOT AND PIVOT INDEX C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IPVT = JMINN ROWMAX = DABS(BOTBLK(JMINN,J)) LOOP = JMINN+1 DO 330 I = LOOP, NRWBOT TEMPIV = DABS(BOTBLK(I,J)) IF (TEMPIV.LE.ROWMAX) GO TO 330 IPVT = I ROWMAX = TEMPIV 330 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C TEST FOR SINGULARITY: C C IF SINGULAR THEN TERMINATE AT 1000; C ELSE CONTINUE. C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF (PIVMAX+ROWMAX.EQ.PIVMAX) GO TO 410 PIVMAX = DMAX1(ROWMAX,PIVMAX) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C IF NECESSARY INTERCHANGE ROWS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INCRJ = INCR+J PIVOT(INCRJ) = INCR+IPVT+NRWTOP IF (IPVT.EQ.JMINN) GO TO 350 DO 340 L = J, NOVRLP SWAP = BOTBLK(IPVT,L) BOTBLK(IPVT,L) = BOTBLK(JMINN,L) BOTBLK(JMINN,L) = SWAP 340 CONTINUE 350 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C COMPUTE MULTIPLIERS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ROWPIV = BOTBLK(JMINN,J) DO 360 I = LOOP, NRWBOT BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV 360 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C PERFORM ROW ELIMINATION WITH COLUMN INDEXING C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 380 L = JPLUS1, NOVRLP ROWMLT = BOTBLK(JMINN,L) DO 370 I = LOOP, NRWBOT BOTBLK(I,L) = BOTBLK(I,L)-ROWMLT*BOTBLK(I,J) 370 CONTINUE 380 CONTINUE 390 CONTINUE 400 CONTINUE C C*************************************************************** C C DONE PROVIDED THE LAST ELEMENT IS NOT ZERO C C*************************************************************** C IF (PIVMAX+DABS(BOTBLK(NRWBOT,NOVRLP)).NE.PIVMAX) RETURN C C*************************************************************** C C **** MATRIX IS SINGULAR - SET IFLAG = - 1. C TERMINATE AT 1000. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C 410 CONTINUE IFLAG = -1 RETURN END SUBROUTINE CRSLVE (TOPBLK,NRWTOP,NOVRLP,ARRAY,NRWBLK,NCLBLK,NBLOKS * ,BOTBLK,NRWBOT,PIVOT,B,X) C C*************************************************************** C C C R S L V E SOLVES THE LINEAR SYSTEM C A*X = B C USING THE DECOMPOSITION ALREADY GENERATED IN C R D C M P. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ***** PARAMETERS ***** C C *** ON ENTRY ... C C TOPBLK - DOUBLE PRECISION(NRWTOP,NOVRLP) C OUTPUT FROM C R D C M P C C NOVRLP - INTEGER C THE NUMBER OF COLUMNS IN WHICH SUCC- C ESSIVE BLOCKS OVERLAP, WHERE C NOVRLP = NRWTOP + NRWBOT C C NRWTOP - INTEGER C NUMBER OF ROWS IN THE BLOCK TOPBLK C C ARRAY - DOUBLE PRECISION(NRWBLK,NCLBLK,NBLOKS) C OUTPUT FROM C R D C M P C C NRWBLK - INTEGER C NUMBER OF ROWS IN K-TH BLOCK C C NCLBLK - INTEGER C NUMBER OF COLUMNS IN K-TH BLOCK C C NBLOKS - INTEGER C NUMBER OF NRWBLK BY NCLBLK BLOCKS IN C THE MATRIX A C C BOTBLK - DOUBLE PRECISION(NRWBOT,NOVRLP) C OUTPUT FROM C R D C M P C C NRWBOT - INTEGER C NUMBER OF ROWS IN THE BLOCK BOTBLK C C PIVOT - INTEGER(N) C THE PIVOT VECTOR FROM C R D C M P C C B - DOUBLE PRECISION(N) C THE RIGHT HAND SIDE VECTOR C C X - DOUBLE PRECISION(N) C WORK SPACE C C *** ON RETURN ... C C X - DOUBLE PRECISION(N) C THE SOLUTION VECTOR C C*************************************************************** C IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER PIVOT(1) DIMENSION TOPBLK(NRWTOP,1),ARRAY(NRWBLK,NCLBLK,1),BOTBLK(NRWBOT,1) * ,B(1),X(1) C C*************************************************************** C C **** DEFINE THE CONSDTANTS USED THROUGHOUT **** C C*************************************************************** C NRWTP1 = NRWTOP+1 NRWBK1 = NRWBLK+1 NVRLP1 = NOVRLP+1 NRWTP0 = NRWTOP-1 NRWBT1 = NRWBOT+1 NROWEL = NRWBLK-NRWTOP NRWEL1 = NROWEL+1 NVRLP0 = NOVRLP-1 NBLKS1 = NBLOKS+1 NBKTOP = NRWBLK+NRWTOP C C*************************************************************** C C **** FORWARD RECURSION **** C C*************************************************************** C C *** FIRST, IN TOPBLK.... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FORWARD SOLUTION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 30 J = 1, NRWTOP X(J) = B(J)/TOPBLK(J,J) IF (J.EQ.NRWTOP) GO TO 20 XJ = -X(J) LOOP = J+1 DO 10 I = LOOP, NRWTOP B(I) = B(I)+TOPBLK(I,J)*XJ 10 CONTINUE 20 CONTINUE 30 CONTINUE C C ******************************************************** C C *** IN EACH BLOCK ARRAY(,,K).... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INCR = 0 DO 120 K = 1, NBLOKS INCRTP = INCR+NRWTOP C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FORWARD MODIFICATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 50 J = 1, NRWTOP INCRJ = INCR+J XINCRJ = -X(INCRJ) DO 40 I = 1, NRWBLK INCRI = INCRTP+I B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ 40 CONTINUE 50 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FORWARD ELIMINATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 80 J = NRWTP1, NRWBLK INCRJ = INCR+J JPIVOT = PIVOT(INCRJ) IF (JPIVOT.EQ.INCRJ) GO TO 60 SWAP = B(INCRJ) B(INCRJ) = B(JPIVOT) B(JPIVOT) = SWAP 60 CONTINUE BINCRJ = -B(INCRJ) LOOP = J-NRWTP0 DO 70 I = LOOP, NRWBLK INCRI = INCRTP+I B(INCRI) = B(INCRI)+ARRAY(I,J,K)*BINCRJ 70 CONTINUE 80 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FORWARD SOLUTION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 110 J = NRWBK1, NBKTOP INCRJ = INCR+J JRWTOP = J-NRWTOP X(INCRJ) = B(INCRJ)/ARRAY(JRWTOP,J,K) IF (J.EQ.NBKTOP) GO TO 100 XINCRJ = -X(INCRJ) LOOP = J-NRWTP0 DO 90 I = LOOP, NRWBLK INCRI = INCRTP+I B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ 90 CONTINUE 100 CONTINUE 110 CONTINUE INCR = INCR+NRWBLK 120 CONTINUE C C ******************************************************** C C *** FINALLY, IN BOTBLK.... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FORWARD MODIFICATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INCRTP = INCR+NRWTOP DO 140 J = 1, NRWTOP INCRJ = INCR+J XINCRJ = -X(INCRJ) DO 130 I = 1, NRWBOT INCRI = INCRTP+I B(INCRI) = B(INCRI)+BOTBLK(I,J)*XINCRJ 130 CONTINUE 140 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C FORWARD ELIMINATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF (NRWBOT.EQ.1) GO TO 180 DO 170 J = NRWTP1, NVRLP0 INCRJ = INCR+J JPIVOT = PIVOT(INCRJ) IF (JPIVOT.EQ.INCRJ) GO TO 150 SWAP = B(INCRJ) B(INCRJ) = B(JPIVOT) B(JPIVOT) = SWAP 150 CONTINUE BINCRJ = -B(INCRJ) LOOP = J-NRWTP0 DO 160 I = LOOP, NRWBOT INCRI = INCRTP+I B(INCRI) = B(INCRI)+BOTBLK(I,J)*BINCRJ 160 CONTINUE 170 CONTINUE 180 CONTINUE C C*************************************************************** C C **** BACKWARD RECURSION **** C C*************************************************************** C C *** FIRST IN BOTBLK.... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C BACKWARD SOLUTION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 210 LL = 1, NRWBOT J = NVRLP1-LL INCRJ = INCR+J NRWBTL = NRWBT1-LL X(INCRJ) = B(INCRJ)/BOTBLK(NRWBTL,J) IF (LL.EQ.NRWBOT) GO TO 200 XINCRJ = -X(INCRJ) LOOP = NRWBOT-LL DO 190 I = 1, LOOP INCRI = INCRTP+I B(INCRI) = B(INCRI)+BOTBLK(I,J)*XINCRJ 190 CONTINUE 200 CONTINUE 210 CONTINUE C C ******************************************************** C C *** THEN IN EACH BLOCK ARRAY(,,K).... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 300 L = 1, NBLOKS C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C BACKWARD ELIMINATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C K = NBLKS1-L INCR = INCR-NRWBLK DO 240 L1 = NRWEL1, NRWBLK I = NRWBLK+NRWEL1-L1 IPLUSN = I+NRWTOP LOOP = IPLUSN+1 INCRN = INCR+IPLUSN DOTPRD = X(INCRN) DO 220 J = LOOP, NCLBLK INCRJ = INCR+J DOTPRD = DOTPRD-ARRAY(I,J,K)*X(INCRJ) 220 CONTINUE X(INCRN) = DOTPRD IPVTN = PIVOT(INCRN) IF (INCRN.EQ.IPVTN) GO TO 230 SWAP = X(INCRN) X(INCRN) = X(IPVTN) X(IPVTN) = SWAP 230 CONTINUE 240 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C BACKWARD MODIFICATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INCRTP = INCR+NRWTOP DO 260 J = NRWBK1, NCLBLK INCRJ = INCR+J XINCRJ = -X(INCRJ) DO 250 I = 1, NROWEL INCRI = INCRTP+I B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ 250 CONTINUE 260 CONTINUE C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C BACKWARD SOLUTION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 290 LL = 1, NROWEL J = NRWBK1-LL INCRJ = INCR+J NRWELL = NRWEL1-LL X(INCRJ) = B(INCRJ)/ARRAY(NRWELL,J,K) IF (LL.EQ.NROWEL) GO TO 280 XINCRJ = -X(INCRJ) LOOP = NROWEL-LL DO 270 I = 1, LOOP INCRI = INCRTP+I B(INCRI) = B(INCRI)+ARRAY(I,J,K)*XINCRJ 270 CONTINUE 280 CONTINUE 290 CONTINUE 300 CONTINUE C C ******************************************************** C C *** IN TOPBLK FINISH WITH.... C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C BACKWARD ELIMINATION C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DO 330 L = 1, NRWTOP I = NRWTP1-L LOOP = I+1 DOTPRD = X(I) DO 310 J = LOOP, NOVRLP DOTPRD = DOTPRD-TOPBLK(I,J)*X(J) 310 CONTINUE X(I) = DOTPRD IPVTI = PIVOT(I) IF (I.EQ.IPVTI) GO TO 320 SWAP = X(I) X(I) = X(IPVTI) X(IPVTI) = SWAP 320 CONTINUE 330 CONTINUE RETURN END subroutine lufac(n, ndim, a, ip, ier) implicit double precision (a-h,o-z) dimension a(ndim,n), ip(n) intrinsic abs * blas: daxpy, dscal, dswap. idamax parameter ( zero = 0.0d+0, one = 1.0d+0 ) * The subroutine lufac is a very simple code to compute the * LU decomposition (with partial pivoting) of the n by n matrix a. * The LU factors are overwritten on a. The integer array ip * reflects the pairwise interchanges performed. note that ip(k) * therefore does not give the index in the original array of * the k-th pivot. * On exit, the error flag ier is zero when no zero pivots are * encountered. Otherwise, ier is equal to the index of the * step at which a zero pivot occurred. ier = 0 ip(n) = 0 * Begin loop over columns 1 through n-1. k is the current * column index. do 100 k = 1, n-1 * Find the row index ipiv of the element of largest magnitude in * column k. ipiv = k-1 + idamax(n-k+1, a(k,k), 1) piv = a(ipiv,k) if (piv .eq. zero) then ier = k return endif ip(k) = ipiv * Perform interchanges if necessary. if (ipiv .ne. k) then call dswap(n-k+1, a(ipiv,k), ndim, a(k,k), ndim) endif * Save the (negative) multipliers in the subdiagonal elements of * column k. call dscal(n-k, (-one/piv), a(k+1,k), 1) * Update the remaining matrix. Note that a(i,k) now contains * the negative multipliers. do 50 j = k+1, n call daxpy(n-k, a(k,j), a(k+1,k), 1, a(k+1,j), 1) 50 continue * End of loop over columns. 100 continue if (a(n,n).eq.zero) ier = n return end subroutine lusol (n, ndim, a, ip, b, x) implicit double precision (a-h, o-z) dimension a(ndim,n), ip(n), b(n), x(n) * blas: daxpy, dcopy * The subroutine lusol is a simple-minded routine to solve a * linear system whose LU factors have been computed by lufac. * On entry, the matrix a should contain the LU factors, and * ip should contain the interchange array constructed by lufac. * Copy the right-hand side b into x. call dcopy (n, b, 1, x, 1) * Forward solution with l (unit lower-triangular factor), which * is stored in the strict lower triangle of a. do 20 k = 1, n-1 ipiv = ip(k) if (ipiv .ne. k) then tem = x(ipiv) x(ipiv) = x(k) x(k) = tem endif call daxpy ( n-k, x(k), a(k+1,k), 1, x(k+1), 1 ) 20 continue * Backward solution with u (upper-triangular factor), which is stored * in the upper triangle of a. do 40 kb = n, 1, -1 x(kb) = x(kb)/a(kb,kb) call daxpy(kb-1, (-x(kb)), a(1,kb), 1, x(1), 1) 40 continue return end subroutine dcopy ( n, x, incx, y, incy ) implicit double precision (a-h,o-z) integer n, incx, incy dimension x( * ), y( * ) c dcopy performs the operation c c y := x c c nag fortran 77 version of the blas routine dcopy . c nag fortran 77 o( n ) basic linear algebra routine. c c -- written on 26-november-1982. c sven hammarling, nag central office. integer i , ix , iy if( n.lt.1 )return if( ( incx.eq.incy ).and.( incy.gt.0 ) )then do 10, iy = 1, 1 + ( n - 1 )*incy, incy y( iy ) = x( iy ) 10 continue else if( incx.ge.0 )then ix = 1 else ix = 1 - ( n - 1 )*incx end if if( incy.gt.0 )then do 20, iy = 1, 1 + ( n - 1 )*incy, incy y( iy ) = x( ix ) ix = ix + incx 20 continue else iy = 1 - ( n - 1 )*incy do 30, i = 1, n y( iy ) = x( ix ) iy = iy + incy ix = ix + incx 30 continue end if end if return * end of dcopy . end subroutine daxpy ( n, alpha, x, incx, y, incy ) implicit double precision (a-h,o-z) integer n, incx, incy dimension x( * ), y( * ) c daxpy performs the operation c c y := alpha*x + y c c c modified nag fortran 77 version of the blas routine daxpy . c c -- written on 3-september-1982. c sven hammarling, nag central office. integer i , ix , iy parameter ( zero = 0.0+0 ) if( n .lt.1 )return if( alpha.eq.zero )return if( ( incx.eq.incy ).and.( incx.gt.0 ) )then do 10, ix = 1, 1 + ( n - 1 )*incx, incx y( ix ) = alpha*x( ix ) + y( ix ) 10 continue else if( incy.ge.0 )then iy = 1 else iy = 1 - ( n - 1 )*incy end if if( incx.gt.0 )then do 20, ix = 1, 1 + ( n - 1 )*incx, incx y( iy ) = alpha*x( ix ) + y( iy ) iy = iy + incy 20 continue else ix = 1 - ( n - 1 )*incx do 30, i = 1, n y( iy ) = alpha*x( ix ) + y( iy ) ix = ix + incx iy = iy + incy 30 continue end if end if return * end of daxpy . end double precision function ddot ( n, x, incx, y, incy ) implicit double precision (a-h,o-z) integer n, incx, incy dimension x( * ), y( * ) c ddot returns the value c c ddot = x'y c c c modified nag fortran 77 version of the blas routine ddot . c c -- written on 21-september-1982. c sven hammarling, nag central office. integer i , ix , iy parameter ( zero = 0.0d+0 ) sum = zero if( n.ge.1 )then if( ( incx.eq.incy ).and.( incx.gt.0 ) )then do 10, ix = 1, 1 + ( n - 1 )*incx, incx sum = sum + x( ix )*y( ix ) 10 continue else if( incy.ge.0 )then iy = 1 else iy = 1 - ( n - 1 )*incy end if if( incx.gt.0 )then do 20, ix = 1, 1 + ( n - 1 )*incx, incx sum = sum + x( ix )*y( iy ) iy = iy + incy 20 continue else ix = 1 - ( n - 1 )*incx do 30, i = 1, n sum = sum + x( ix )*y( iy ) ix = ix + incx iy = iy + incy 30 continue end if end if end if ddot = sum return * end of ddot . end subroutine dscal ( n, alpha, x, incx ) implicit double precision (a-h,o-z) integer n, incx dimension x( * ) c dscal performs the operation c c x := alpha*x c c c modified nag fortran 77 version of the blas routine dscal . c c -- written on 26-november-1982. c sven hammarling, nag central office. integer ix parameter ( one = 1.0d+0, zero = 0.0d+0 ) if( n.ge.1 )then if( alpha.eq.zero )then do 10, ix = 1, 1 + ( n - 1 )*incx, incx x( ix ) = zero 10 continue else if( alpha.eq.( -one ) )then do 20, ix = 1, 1 + ( n - 1 )*incx, incx x( ix ) = -x( ix ) 20 continue else if( alpha.ne.one )then do 30, ix = 1, 1 + ( n - 1 )*incx, incx x( ix ) = alpha*x( ix ) 30 continue end if end if return * end of dscal . end subroutine dswap ( n, x, incx, y, incy ) implicit double precision (a-h,o-z) integer n, incx, incy dimension x( * ), y( * ) c dswap performs the operations c c temp := x, x := y, y := temp. c c c modified nag fortran 77 version of the blas routine dswap . c c -- written on 26-november-1982. c sven hammarling, nag central office. integer i , ix , iy if( n.lt.1 )return if( ( incx.eq.incy ).and.( incy.gt.0 ) )then do 10, iy = 1, 1 + ( n - 1 )*incy, incy temp = x( iy ) x( iy ) = y( iy ) y( iy ) = temp 10 continue else if( incx.ge.0 )then ix = 1 else ix = 1 - ( n - 1 )*incx end if if( incy.gt.0 )then do 20, iy = 1, 1 + ( n - 1 )*incy, incy temp = x( ix ) x( ix ) = y( iy ) y( iy ) = temp ix = ix + incx 20 continue else iy = 1 - ( n - 1 )*incy do 30, i = 1, n temp = x( ix ) x( ix ) = y( iy ) y( iy ) = temp iy = iy + incy ix = ix + incx 30 continue end if end if return * end of dswap . end integer function idamax( n, x, incx ) implicit double precision (a-h,o-z) integer n, incx dimension x( * ) c idamax returns the smallest value of i such that c c abs( x( i ) ) = max( abs( x( j ) ) ) c j c c nag fortran 77 version of the blas routine idamax. c nag fortran 77 o( n ) basic linear algebra routine. c c -- written on 31-may-1983. c sven hammarling, nag central office. intrinsic abs integer i , imax , ix if( n.lt.1 )then idamax = 0 return end if imax = 1 if( n.gt.1 )then xmax = abs( x( 1 ) ) ix = 1 do 10, i = 2, n ix = ix + incx if( xmax.lt.abs( x( ix ) ) )then xmax = abs( x( ix ) ) imax = i end if 10 continue end if idamax = imax return * end of idamax. end subroutine dload ( n, const, x, incx ) implicit double precision (a-h,o-z) dimension x(*) c c dload performs the operation c c x = const*e, e' = ( 1 1 ... 1 ). c c c nag fortran 77 o( n ) basic linear algebra routine. c c -- written on 22-september-1983. c sven hammarling, nag central office. c parameter ( zero = 0.0d+0 ) if( n.lt.1 )return if( const.ne.zero )then do 10, ix = 1, 1 + ( n - 1 )*incx, incx x( ix ) = const 10 continue else do 20, ix = 1, 1 + ( n - 1 )*incx, incx x( ix ) = zero 20 continue end if return * end of dload . end subroutine maxpy ( nrow, ncol, alpha, xmat, nrowy, ymat ) implicit double precision (a-h,o-z) dimension xmat(nrow, ncol), ymat(nrowy, ncol) * Subroutine maxpy takes as input the scalar alpha and two matrices, * xmat and ymat. xmat has declared row dimension nrow, and * ymat has declared row dimension nrowy, but both are * conceptually nrow by ncol. * On output, (new ymat) is alpha*xmat+ (old ymat), by analogy * with the vector blas routine saxpy. do 100 j = 1, ncol do 100 i = 1, nrow ymat(i,j) = ymat(i,j) + alpha*xmat(i,j) 100 continue return end subroutine matcop( nrow1, nrow2, nrow, ncol, xmat1, xmat2 ) implicit double precision (a-h,o-z) dimension xmat1(nrow1, ncol), xmat2(nrow2, ncol) * Given 2 matrices xmat1 and xmat2, where xmat1 has declared * row dimension nrow1, xmat2 has declared row dimension nrow2, * and both have column dimension ncol, the routine matcop copies * rows 1 through nrow, and columns 1 through ncol from xmat1 into * xmat2. if (nrow .le. 0 .or. ncol .le. 0) return do 100 j = 1, ncol do 100 i = 1, nrow xmat2(i,j) = xmat1(i,j) 100 continue return end subroutine mtload( nrow, ncol, const, nrowx, xmat ) implicit double precision (a-h,o-z) dimension xmat(nrowx, ncol) * mtload sets elements 1 through nrow, 1 through ncol, of the * matrix xmat (whose declared row dimension is nrowx) to the * scalar value const. if (nrow .le. 0 .or. ncol .le. 0) return do 100 j = 1, ncol do 100 i = 1, nrow xmat(i,j) = const 100 continue return end subroutine mssq ( nrow, ncol, xmat, scale, sumsq ) implicit double precision (a-h,o-z) dimension xmat(nrow, *) * Given the nrow by ncol matrix xmat, mssq returns values * scale and sumsq such that * (scale**2) * sumsq = sum of squares of xmat(i,j), * where scale = max abs(xmat(i,j)). * mssq is a stripped-down matrix version of the blas routine sssq. intrinsic abs parameter ( one = 1.0d+0, zero = 0.0d+0 ) scale = zero sumsq = one if( nrow.ge.1 .and. ncol .ge. 1) then do 10 i = 1, nrow do 10 j = 1, ncol if( xmat(i,j) .ne. zero )then absxij = abs(xmat(i,j)) if( scale .lt. absxij ) then sumsq = one + sumsq* (scale/absxij)**2 scale = absxij else sumsq = sumsq + (absxij/scale)**2 end if end if 10 continue end if return end subroutine dssq ( n, x, incx, scale, sumsq ) implicit double precision (a-h,o-z) integer n, incx dimension x( * ) * Given the n-vector x, dssq returns values scale and sumsq such that * (scale**2) * sumsq = sum of squares of x(i), * where scale = max abs(x(i)). intrinsic abs parameter ( one = 1.0d+0, zero = 0.0d+0 ) scale = zero sumsq = one if( n.ge.1 )then do 10, ix = 1, 1 + ( n - 1 )*incx, incx if( x( ix ).ne.zero )then absxi = abs( x( ix ) ) if( scale.lt.absxi )then sumsq = one + sumsq*( scale/absxi )**2 scale = absxi else sumsq = sumsq + ( absxi/scale )**2 end if end if 10 continue end if return * end of dssq . end subroutine bvpsol(ncomp, nmsh, nlbc, aleft, aright, * nfxpnt, fixpnt, * ntol, ltol, tol, nmax, linear, giveu, givmsh, * xx, nudim, u, defexp, defimp, def, delu, rhs, fval, * topblk, botblk, ajac, bhold, chold, ipvblk, ipivlu, * uint, ftmp, tmprhs, dftmp1, dftmp2, dgtm, * utrial, rhstri, xmerit, xxold, uold, usave, * tmp, dsq, dexr, ratdc, rerr, * etest6, etest8, ermx, ihcomp, irefin, * def6, def8, fsub, dfsub, gsub, dgsub, iflbvp) implicit double precision (a-h,o-z) dimension fixpnt(*), ltol(ntol), tol(ntol) dimension xx(*), u(nudim, *) dimension defexp(ncomp,*), defimp(ncomp,*), def(ncomp,*) dimension delu(ncomp, *), rhs(*), fval(ncomp,*) dimension topblk(nlbc,*), botblk(ncomp-nlbc,*) dimension ajac(ncomp, 2*ncomp, *) dimension bhold(ncomp, ncomp, *), chold(ncomp, ncomp, *) dimension ipivlu(*), ipvblk(*) dimension uint(ncomp), ftmp(ncomp) dimension dgtm(ncomp), tmprhs(*) dimension dftmp1(ncomp, ncomp), dftmp2(ncomp, ncomp) dimension utrial(ncomp,*), rhstri(*) dimension xmerit(ncomp, *) dimension xxold(*), uold(ncomp,*), usave(ncomp,*) dimension tmp(ncomp,8) dimension dsq(ncomp,ncomp), dexr(ncomp) dimension ratdc(*), rerr(ncomp,*) dimension etest6(*), etest8(*), ermx(*) dimension ihcomp(*), irefin(*) dimension def6(ncomp,*), def8(ncomp,*) logical linear, giveu, givmsh, double external fsub, dfsub, gsub, dgsub common/mchprs/flmin, flmax, epsmch logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 intrinsic max logical smooth, succes, strctr, trst6, reaft6 logical onto6, onto8, ludone, rhsgiv, maxmsh logical first4, first8 logical mchset save mchset parameter (zero = 0.0d+0, one = 1.0d+0) parameter (third = 0.33d+0, fourth = 0.25d+0) parameter (quan6 = 0.1d+0 ) * blas: dload * double precision d1mach data mchset/.true./ data fxfct/10.0d+0/ data maxmsh/.false./ if (mchset) then flmin = d1mach(1) flmax = d1mach(2) epsmch = d1mach(3) if (pdebug) write(6,901) epsmch mchset = .false. endif * The routine stcons calculates integration constants stored in * labeled common consts. call stcons * Set up arrays for the error tests. if (.not. linear) then call dload(ntol, one, etest6, 1) else do 10 i = 1, ntol etest6(i) = one/max(quan6, tol(i)**third) 10 continue endif nmold = 1 smooth = .false. strctr = .false. trst6 = .true. reaft6 = .false. numbig = 0 nummed = 0 first4 = .true. first8 = .true. * * If givmsh is .true., the initial number of mesh points must be * provided by the user in nmsh, and the mesh points must be * contained in the array xx (of dimension nmsh). * Otherwise, nmsh is set to its default value, and a * uniform initial mesh is created. if (.not. giveu .and. .not. givmsh) then nmsh = nminit if (nmsh .lt. nfxpnt+2) nmsh = nfxpnt + 2 call unimsh(nmsh, aleft, aright, nfxpnt, fixpnt, xx) endif if (pdebug) then write(6,902) call sprt(nmsh, xx) endif if (.not. giveu) call initu(ncomp, nmsh, xx, nudim, u) ***** top of logic for 4th order solution **** 400 continue if (iprint .eq. 1) write(6,903) nmsh * Set the def (deferred correction) array to zero. call mtload(ncomp, nmsh-1, zero, ncomp, def) iorder = 4 * The routine fneval calls fsub at the mesh xx and the * solution u, and saves the values in the array fval. call fneval(ncomp, nmsh, xx, nudim, u, fval, fsub) * Try to compute a 4th order solution by solving a system of nonlinear * equations. if (linear) then ludone = .false. call lineq( ncomp, nmsh, nlbc, * ludone, xx, nudim, u, def, * delu, rhs, fval, * uint, ftmp, dftmp1, dftmp2, dgtm, tmprhs, * ajac, topblk, botblk, bhold, chold, ipvblk, * fsub, dfsub, gsub, dgsub, iflnwt) * Call fneval to evaluate the fval array at the new solution u. * (Such a call is not necessary for the nonlinear case because * fval is called within newteq for the new u.) call fneval(ncomp, nmsh, xx, nudim, u, fval, fsub) else rhsgiv = .false. call newteq(ncomp, nmsh, nlbc, * rhsgiv, ntol, ltol, tol, * xx, nudim, u, def, * delu, rhs, fval, * utrial, rhstri, * uint, ftmp, dftmp1, dftmp2, dgtm, tmprhs, xmerit, * ajac, topblk, botblk, bhold, chold, ipvblk, * fsub, dfsub, gsub, dgsub, itnwt, iflnwt) endif if (iflnwt .eq. 0) then call conv4( ncomp, nmsh, ntol, ltol, tol, linear, nmax, * xx, nudim, u, defexp, defimp, def, fval, * tmp, bhold, chold, dsq, dexr, usave, * ratdc, rerr, ipivlu, nmold, xxold, * smooth, reaft6, onto6, strctr, trst6, double, * fsub, maxmsh, succes, first4) if (pdebug .and. .not. onto6) write (6,904) else succes = .false. onto6 = .false. call fail4( ncomp, nmsh, nlbc, ntol, ltol, * xx, nudim, u, rhs, linear, nmax, * nmold, xxold, uold, ratdc, * iorder, iflnwt, itnwt, double, maxmsh, * numbig, nummed) endif if (succes) then iflbvp = 0 return elseif (maxmsh) then go to 900 elseif (.not. onto6) then go to 400 endif * To reach here, onto6 must be .true. **** logic for 6th order **** if (iprint .eq. 1) write(6,905) * Save the 4th order solution on this mesh in uold. call matcop(nudim, ncomp, ncomp, nmsh, u, uold) iorder = 6 if (linear) then call lineq( ncomp, nmsh, nlbc, * ludone, xx, nudim, u, def, * delu, rhs, fval, * uint, ftmp, dftmp1, dftmp2, dgtm, tmprhs, * ajac, topblk, botblk, chold, bhold, ipvblk, * fsub, dfsub, gsub, dgsub, iflnwt) else call fixjac(ncomp, nmsh, nlbc, * iorder, ntol, ltol, tol, * xx, nudim, u, def, def, delu, * rhs, fval, utrial, rhstri, * rnsq, uint, ftmp, tmprhs, * ajac, topblk, botblk, ipvblk, * fsub, gsub, iflnwt) * If the fixed Jacobian iterations fail but rnsq is small, * try a Newton procedure. Set rhsgiv to indicate that * the right-hand side and fval have already been evaluated * at the given u. if (iflnwt .eq. -3 .and. rnsq .lt. fxfct*epsmch) then rhsgiv = .true. call newteq(ncomp, nmsh, nlbc, * rhsgiv, ntol, ltol, tol, * xx, nudim, u, def, * delu, rhs, fval, * utrial, rhstri, * uint, ftmp, dftmp1, dftmp2, dgtm, tmprhs, xmerit, * ajac, topblk, botblk, bhold, chold, ipvblk, * fsub, dfsub, gsub, dgsub, iter, iflnwt) endif endif if (iflnwt .eq. 0) then call conv6(ncomp, nmsh, ntol, ltol, tol, * nudim, u, uold, etest6, err6, * trst6, onto8, reaft6, succes) else onto8 = .false. call fail6( ncomp, nmsh, nlbc, ntol, ltol, tol, * nfxpnt, fixpnt, iorder, nmax, * xx, nudim, u, rhs, usave, xxold, uold, nmold, * ihcomp, irefin, * rerr, ermx, ratdc, * reaft6, double, succes, maxmsh, * numbig, nummed) endif if (succes) then iflbvp = 0 return elseif (maxmsh) then go to 900 elseif (.not. onto8) then go to 400 endif ***** logic for trying to calculate 8th order solution ***** if (iprint .eq. 1) write(6,906) call matcop(nudim, ncomp, ncomp, nmsh, u, uold) * Save the old deferred correction vector def in def6. call matcop(ncomp, ncomp, ncomp, nmsh-1, def, def6) * For linear problems, calculate the fval array for the * new solution u. if (linear) call fneval(ncomp, nmsh, xx, nudim, u, fval, fsub) * Calculate 8th order deferred corrections (the array def8). call df8cal (ncomp, nmsh, xx, nudim, u, fval, def8, * tmp, fsub) * For linear problems, the def array is the def8 array. * For nonlinear problems, add the def8 array to the * already-calculated def array. if (linear) then call matcop(ncomp, ncomp, ncomp, nmsh-1, def8, def) else call maxpy(ncomp, nmsh-1, one, def8, ncomp, def) endif iorder = 8 if (linear) then call lineq( ncomp, nmsh, nlbc, * ludone, xx, nudim, u, def, * delu, rhs, fval, * uint, ftmp, dftmp1, dftmp2, dgtm, tmprhs, * ajac, topblk, botblk, chold, bhold, ipvblk, * fsub, dfsub, gsub, dgsub, iflnwt) else call fixjac(ncomp, nmsh, nlbc, * iorder, ntol, ltol, tol, * xx, nudim, u, def, def8, delu, * rhs, fval, utrial, rhstri, * rnsq, uint, ftmp, tmprhs, * ajac, topblk, botblk, ipvblk, * fsub, gsub, iflnwt) * If the fixed Jacobian iterations fail but rnsq is small, * try a Newton procedure. Set rhsgiv to indicate that * the right-hand side and fval have already been evaluated * at the given u. if (iflnwt .eq. -3 .and. rnsq .lt. fxfct*epsmch) then rhsgiv = .true. call newteq(ncomp, nmsh, nlbc, * rhsgiv, ntol, ltol, tol, * xx, nudim, u, def, * delu, rhs, fval, * utrial, rhstri, * uint, ftmp, dftmp1, dftmp2, dgtm, tmprhs, xmerit, * ajac, topblk, botblk, bhold, chold, ipvblk, * fsub, dfsub, gsub, dgsub, iter, iflnwt) endif endif if (iflnwt .eq. 0) then call conv8( ncomp, nmsh, ntol, ltol, tol, * nfxpnt, fixpnt, linear, nmax, * xx, nudim, u, def, def6, def8, uold, * ihcomp, irefin, ermx, err6, * etest8, strctr, * double, nmold, xxold, maxmsh, succes, first8) else succes = .false. call fail8(ncomp, nmsh, nfxpnt, fixpnt, nmax, * ntol, ltol, tol, nmold, * xx, nudim, u, def6, xxold, uold, * ihcomp, irefin, ermx, double, maxmsh) endif if (maxmsh) then go to 900 elseif (.not. succes) then go to 400 endif * Successful termination. iflbvp = 0 return 900 continue * Error exit---too many mesh points. iflbvp = 1 return 901 format(1h ,'epsmch',1pe10.3) 902 format(1h ,'initial mesh') 903 format(1h ,'start 4th order, nmsh',i5) 904 format(1h ,'do not go on to 6th') 905 format(1h ,'start 6th order') 906 format(1h ,'start 8th order') end subroutine initu(ncomp, nmsh, xx, nudim, u) implicit double precision (a-h,o-z) dimension xx(*), u(nudim, *) logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 * This routine must be provided to reset u after re-meshing * for linear problems or for nonlinear problems * when interpolation of the old solution is not used. * This version sets all elements of u to the constant uval0. if (iprint .ne. -1) write(6,99) uval0 99 format(1h ,'initu, uval0',1pd15.5) call mtload(ncomp, nmsh, uval0, nudim, u) return end block data * This block data routine initializes nminit (the initial number * of mesh points), pdebug (a logical variable indicating whether * debug printout is desired), and uval0 (the initial value for the trial * solution) to their default values. double precision uval0 logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 data nminit/7/ data pdebug/.false./ data iprint/0/ data uval0/0.0d+0/ end subroutine conv4( ncomp, nmsh, ntol, ltol, tol, linear, nmax, * xx, nudim, u, defexp, defimp, def, fval, * tmp, bhold, chold, dsq, dexr, usave, * ratdc, rerr, ipivot, nmold, xxold, * smooth, reaft6, onto6, strctr, trst6, double, * fsub, maxmsh, succes, first4) implicit double precision (a-h,o-z) dimension ltol(ntol), ipivot(*) dimension xx(*), u(nudim,*), tol(ntol) dimension defexp(ncomp,*), defimp(ncomp,*), def(ncomp,*) dimension fval(ncomp,*), tmp(ncomp,4) dimension bhold(ncomp,ncomp,*), chold(ncomp, ncomp,*) dimension dsq(ncomp,ncomp), dexr(ncomp), usave(ncomp,*) dimension xxold(*) dimension ratdc(*), rerr(ncomp,*) logical linear, smooth, reaft6, onto6, strctr, trst6 logical double, succes, maxmsh, first4 external fsub logical callrt, oscchk, reposs, savedu logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 common/mchprs/flmin, flmax, epsmch parameter (hundth = 1.0d+5, rerfct = 1.5d+0) parameter (power = 1.0d+0/6.0d+0, one = 1.0d+0) parameter (zero = 0.0d+0, huge = 1.0d+30) intrinsic abs, max save dfold, oldrt1, savedu, reposs logical adjrer * The Newton iteration converged for a 4th order solution. if (iprint .eq. 1) write(6,901) if (first4) then dfold = zero oldrt1 = huge savedu = .false. reposs = .false. first4 = .false. endif succes = .false. maxmsh = .false. * Compute the explicit deferred correction array defexp. * The parameter fval is an ncomp by nmsh array, calculated * by a previous call to fneval for this mesh and solution. call dfexcl( ncomp, nmsh, xx, nudim, u, defexp, fval, * tmp, fsub ) call matcop(ncomp, ncomp, ncomp, nmsh-1, defexp, def) * If the problem has been declared as smooth, or we might wish * to perform Richardson extrapolation after trying to calculate a * 6th order solution, just go on to the 6th order logic (don't * bother to calculate implicit deferred corrections). if ( smooth .or. reaft6) then if (smooth .and. pdebug) write(6,902) if (reaft6 .and. pdebug) write(6,903) onto6 = .true. return endif * Compute the cheap implicit deferred correction array defimp. * The array chold must be unchanged from a call to jaccal * for this mesh and solution. * The temporary arrays dsq and dexr are calculated inside dfimcl. call dfimcl( ncomp, nmsh, defexp, chold, dsq, dexr, * ipivot, defimp) * Call dccal to calculate: dfexmx, the maximum-magnitude element * of defexp in components for which tolerances are specified; * incmp, inmsh, and intol, the indices of the component, * mesh interval, and tolerance of dfexmx; derivm, the * maximum-magnitude element of fval(incmp,*) for all mesh * points; dfimmx, the maximum-magnitude element of defimp in * component incmp; the ratios rat1 and rat2; and the array * ratdc (used in osc). dfctol = hundth*epsmch call dccal( ncomp, nmsh, ntol, ltol, * defexp, defimp, dfctol, fval, * ratdc, dfexmx, incmp, inmsh, intol, * derivm, dfimmx, rat1, rat2) * decid4 sets logical flags that determine the next step of the * algorithm. call decid4(linear, rat1, rat2, dfexmx, dfimmx, * derivm, dfold, tol(intol), oldrt1, * onto6, smooth, callrt, strctr, oscchk, double, reposs) if (pdebug) then if (smooth) write(6,904) if (callrt) write(6,905) if (oscchk) write(6,906) if (strctr) write(6,907) if (reposs) write(6,908) if (savedu) write(6,909) if (onto6) write(6,910) write(6,911) rat1, oldrt1 endif oldrt1 = rat1 dfold = dfexmx if (callrt) then * ratcor calculates a more expensive implicit deferred correction. * The matrix bhold must be unchanged from the last call to jaccal. * If callrt is true, onto6 is always true also. call ratcor( ncomp, nmsh, xx, defimp, bhold, def) elseif (linear) then if (oscchk) then call osc (ncomp, nmsh, dfexmx, incmp, * defexp, ratdc, double, inmsh, onto6, trst6, * smooth) elseif (reposs) then * If reposs (`Richardson extrapolation possible') is true * for two successive 4th order solutions, a special termination * test may be used based on Richardson extrapolation. * If reposs is true for the first time, savedu will * be false; savedu can be true only when reposs is true for a * second consecutive mesh (a doubled version of the first). if (savedu) then * The flag savedu is .true. when the immediately preceding * converged 4th order solution was saved in the array usave, * and the mesh was then doubled. * In this case, the routine rerrvl is called to compute a * Richardson extrapolation (RE) error estimate remax. * The rerr array does not need to be adjusted, so adjrer is false. adjrer = .false. call rerrvl( ncomp, nmsh, nudim, u, usave, ntol, * ltol, rerr, remax, itlmx, adjrer ) if (remax .lt. rerfct*tol(itlmx)) then succes = .true. return endif * end of logic for savedu = .true. endif * Here, reposs is .true., but either we hadn't saved the previous * solution, or else the RE error estimate is not small. * Save u in usave, and set savedu to .true. * Set double to .true. to indicate that the mesh should be doubled. * Set .onto6. to .false. to return to the top of the 4th order * iteration loop. call matcop(nudim, ncomp, ncomp, nmsh, u, usave) savedu = .true. double = .true. onto6 = .false. * end of logic for reposs = .true. endif * end of logic for linear endif if (pdebug .and. reposs .and. .not. onto6) write(6,912) if (.not. onto6) then * NB: onto6 can be false only for linear problems if (double) then call dblmsh (nmsh, nmax, xx, nmold, xxold, maxmsh) else * Refine the mesh near interval inmsh; the number of points * to be added depends on the relative size of dfexmx compared * to u and tol. drat = dfexmx/ * (max(one, abs(u(incmp,inmsh)))*tol(intol)) if (pdebug) write(6,913) drat, u(incmp,inmsh), tol(intol) numadd = drat**power call smpmsh (nmsh, nmax, xx, inmsh, numadd, * nmold, xxold, maxmsh) endif if (.not. maxmsh) call initu(ncomp, nmsh, xx, nudim, u) * end of logic for .not. onto6 endif return 901 format(1h ,'conv4') 902 format(1h ,'smooth') 903 format(1h ,'reaft6') 904 format(1h ,'smooth') 905 format(1h ,'callrt') 906 format(1h ,'oscchk') 907 format(1h ,'strctr') 908 format(1h ,'reposs') 909 format(1h ,'savedu') 910 format(1h ,'onto6 after decid4') 911 format(1h ,'rat1,oldrt1',2(1pe11.3)) 912 format(1h ,'reposs and not onto6') 913 format(1h ,'drat,u,tol',3(1pe11.3)) end subroutine fail4( ncomp, nmsh, nlbc, ntol, ltol, * xx, nudim, u, rhs, linear, nmax, * nmold, xxold, uold, tmwork, * iorder, iflnwt, itnwt, double, maxmsh, * numbig, nummed) implicit double precision (a-h,o-z) dimension ltol(ntol) dimension xx(*), u(nudim, *), rhs(*) dimension xxold(*), uold(ncomp, *), tmwork(*) logical linear, double, maxmsh logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 * The Newton procedure failed to obtain a 4th order solution. if (iprint .eq. 1) write(6,901) maxmsh = .false. if (iflnwt .eq. -1) then * iflnwt = -1 means that the Jacobian was considered singular. * (This is the only possible failure for a linear problem.) call dblmsh (nmsh, nmax, xx, nmold, xxold, maxmsh) call initu(ncomp, nmsh, xx, nudim, u) else * The routine mshref decides how to refine the mesh and then * performs the refinement, either by doubling or based on * the rhs vector at the best point. call mshref(ncomp, nmsh, nlbc, ntol, ltol, * iorder, rhs, tmwork, * nmax, xx, nmold, xxold, double, maxmsh, * numbig, nummed) if (.not. maxmsh) then if (linear .or. itnwt .eq. 0) then call initu(ncomp, nmsh, xx, nudim, u) else * Interpolate the partially converged solution. call matcop(nudim, ncomp, ncomp, nmold, u, uold) call interp(ncomp, nmsh, xx, nudim, u, nmold, * xxold, uold) endif endif * End of logic for failure because of some reason other than * singularity. endif return 901 format(1h ,'fail4') end subroutine conv6(ncomp, nmsh, ntol, ltol, tol, * nudim, u, uold, etest6, err6, * trst6, onto8, reaft6, succes) implicit double precision (a-h,o-z) dimension ltol(ntol) dimension u(nudim,*), tol(ntol) dimension uold(ncomp,*), etest6(*) logical trst6, onto8, reaft6, succes logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 logical errok * The Newton iteration converged for a 6th-order solution. if (iprint .eq. 1) write(6,901) succes = .false. * The logical flag reaft6 is true only when the 6th order solution * failed. Since the 6th order solution converged, reaft6 is false. reaft6 = .false. onto8 = .true. * Calculate the error estimates for the 4th order solution. call errest (ncomp, nmsh, ntol, ltol, tol, * nudim, u, uold, etest6, err6, errok) if (trst6 .and. errok) then succes = .true. return endif return 901 format(1h ,'conv6') end subroutine fail6( ncomp, nmsh, nlbc, ntol, ltol, tol, * nfxpnt, fixpnt, iorder, nmax, * xx, nudim, u, rhs, usave, xxold, uold, nmold, * ihcomp, irefin, rerr, ermx, tmwork, * reaft6, double, succes, maxmsh, * numbig, nummed) implicit double precision (a-h,o-z) dimension ltol(ntol) dimension fixpnt(*), tol(*), rhs(*) dimension xx(*), u(nudim,*), xxold(*) dimension uold(ncomp,*), usave(ncomp,*) dimension ihcomp(*), irefin(*) dimension rerr(ncomp,*), ermx(*), tmwork(*) logical reaft6, double, succes, maxmsh logical adjrer logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 * blas: dcopy parameter (eight = 8.0d+0) * nmpt is the standard number of mesh points to be added to selected * intervals when the mesh is altered based on the distribution * in size of the rhs vector. parameter ( nmpt = 15) * Non-convergence of 6th order. if (iprint .eq. 1) write(6,901) succes = .false. maxmsh = .false. * NB: the problem must be nonlinear. Linear problems will either * fail to converge for 4th order, or else, once they've converged * for 4th order, must converge for 6th order. * Restore the u array to the previous 4th order solution. call matcop(ncomp, nudim, ncomp, nmsh, uold, u) if (reaft6) write(6,9999) 9999 format(1h ,'in fail6, reaft6is true') if (.not.reaft6) write(6,9998) 9998 format(1h ,'in fail6, not reaft6') if (double) write(6,9997) 9997 format(1h ,'in fail6, double is true') if (.not.double) write(6,9996) 9996 format(1h ,'in fail6, not double') if (.not. reaft6 .or. .not. double) then * Here, either * (1) the mesh for which this 6th order solution failed * is not a doubled version of the immediately preceding mesh, or * (2) for the immediately preceding mesh, it is not true * that the 4th order converged and the 6th order failed. * Setting reaft6 to .true. signals that Richardson extrapolation * may be possible if the next 6th order solution fails. When * reaft6 is true, the routine conv4 immediately sets onto6 to true. reaft6 = .true. * Save the current 4th order solution in usave. call matcop(nudim, ncomp, ncomp, nmsh, u, usave) * Use the distribution of sizes in the rhs vector to refine the mesh. call mshref(ncomp, nmsh, nlbc, ntol, ltol, * iorder, rhs, tmwork, * nmax, xx, nmold, xxold, double, maxmsh, * numbig, nummed) if (.not. maxmsh) call interp(ncomp, nmsh, xx, nudim, u, * nmold, xxold, uold) else * Here, reaft6 and double are both true. So for two consecutive * meshes, the 4th order converged and the 6th order failed, * and the second mesh is the double of the first. * Calculate an error estimate from Richardson extrapolation * with the current and previous 4th order solutions. * (usave is the 4th order solution saved from the previous (halved) * mesh.) * Set addrer to .true. to signal that the rerr array should * be adjusted. adjrer = .true. call rerrvl( ncomp, nmsh, nudim, u, usave, ntol, ltol, * rerr, remax, itlmx, adjrer ) write(6,9994) 9994 format(1h ,'***in fail6') write(6,9993) remax, eight*tol(itlmx) 9993 format(1h ,'remax',1pe14.4,5x,'8*tol',1pe14.4) if (remax .lt. eight*tol(itlmx)) then succes = .true. else * Richardson extrapolation did not give sufficient accuracy. * Perform selective mesh refinement on the OLD (NB: old!) mesh * and the old (saved) solution, using the error estimate from * Richardson extrapolation to guide where the mesh points are placed. nmsh = 1 + (nmsh-1)/2 call dcopy(nmsh, xxold, 1, xx, 1) ipow = 4 * The rerr array is overwritten by selmsh. call selmsh(ncomp, nmsh, ntol, ltol, tol, * nfxpnt, fixpnt, ipow, nmax, * xx, ncomp, usave, rerr, irefin, ihcomp, * nmold, xxold, ermx, double, maxmsh) * If double is false on exit from selmsh, the call to selmsh has * produced a different (non-doubled) mesh. Interpolate the * saved solution (from the old mesh) onto the mesh newly created * by selmsh. * NB: Because double is false, we won't try Richardson extrapolation * if the next 4th order converges and 6th order fails. if (.not. maxmsh) then if (.not. double) then call interp(ncomp, nmsh, xx, nudim, u, nmold, * xxold, usave) else * Selective mesh refinement based on the old mesh simply * produced the same mesh that we started with. So now refine * starting with the doubled mesh (from before) and the solution. reaft6 = .true. * Save the solution in usave in case we can carry out Richardson * extrapolation in the same circumstances next time. call matcop(nudim, ncomp, ncomp, nmsh, u, usave) * Use the distribution of sizes in the rhs vector to refine the mesh. call mshref(ncomp, nmsh, nlbc, ntol, ltol, * iorder, rhs, tmwork, * nmax, xx, nmold, xxold, double, maxmsh, * numbig, nummed) if (.not. maxmsh) * call interp(ncomp, nmsh, xx, nudim, u, * nmold, xxold, usave) * end of logic for needing to refine (again) based on the * current mesh endif * end of logic for not too many mesh points endif * end of logic for failure of Richardson extrapolation * to produce a converged solution endif * end of logic for both reaft6 and double being true endif return 901 format(1h ,'fail6') end subroutine conv8( ncomp, nmsh, ntol, ltol, tol, * nfxpnt, fixpnt, linear, nmax, * xx, nudim, u, def, def6, def8, uold, * ihcomp, irefin, ermx, err6, * etest8, strctr, * double, nmold, xxold, maxmsh, succes, first8) implicit double precision (a-h,o-z) dimension ltol(ntol), tol(ntol) dimension fixpnt(*) dimension etest8(ntol) dimension xx(*), u(nudim,*), def(ncomp,*) dimension def6(ncomp,*), def8(ncomp,*), uold(ncomp,*) dimension ihcomp(*), irefin(*) dimension ermx(*), xxold(*) logical linear, strctr, double, maxmsh, succes, first8 logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 intrinsic max logical errok parameter (one = 1.0d+0, fourth = 0.25d+0, quan8 = 0.025d+0) parameter ( efact = 100.0d+0, huge = 1.0d+30 ) * blas: dload save er6old, er8old * The Newton iteration converged for the 8th order solution. if (iprint .eq. 1) write(6,901) if (first8) then er6old = huge er8old = huge first8 = .false. endif if (.not. linear) then call dload(ntol, one, etest8, 1) else do 10 i = 1, ntol etest8(i) = one/max(quan8, tol(i)**fourth) 10 continue endif succes = .false. maxmsh = .false. * Check estimated error. For a nonlinear problem, all components * of etest8 (the ratios used in testing the error) are set to one. * For a linear problem, the components of etest8 are in general * larger than one. But if strctr is .true. and the number of mesh * points decreased, we set the elements of etest8 to one (which * makes a stricter test). if (linear .and. strctr .and. nmsh .lt. nmold) * call dload(ntol, one, etest8, 1) call errest (ncomp, nmsh, ntol, ltol, tol, * nudim, u, uold, etest8, err8, errok) if (errok) then succes = .true. return endif * At this point, the 8th order solution converged, but did not * satisfy the test for termination. if (pdebug) write(6,902) err6, err8, er6old, er8old if (nmsh .lt. nmold. and. * err6 .gt. efact*er6old .and. * err8 .gt. efact*er8old) then * If the number of mesh points decreased and the errors in the * 6th and 8th order solutions did not decrease sufficiently compared * to the previous mesh, double the mesh and go back to try to * calculate a 4th order solution. call dblmsh (nmsh, nmax, xx, nmold, xxold, maxmsh) if (.not. maxmsh) then er6old = err6 er8old = err8 call initu(ncomp, nmsh, xx, nudim, u) endif return endif * Here, we know that * (1) the number of mesh points exceeds that for the previous mesh; or * (2) the number of mesh points decreased, the 6th and 8th order * errors did not satisfy the termination tests, but they did not * increase so much that the mesh needed to be doubled. er6old = err6 er8old = err8 if (err8 .le. err6) then * Perform selective mesh refinement based on the 8th order deferred * corrections. The value of ipow indicates that the error estimate * is of order 6. Then, for a nonlinear problem, interpolate the * latest solution onto the new mesh. ipow = 6 * NB: The array def8 will be overwritten by selmsh. call selmsh(ncomp, nmsh, ntol, ltol, tol, * nfxpnt, fixpnt, ipow, nmax, * xx, nudim, u, def8, irefin, ihcomp, * nmold, xxold, ermx, double, maxmsh) if (.not. maxmsh) then if (linear) then call initu(ncomp, nmsh, xx, nudim, u) else call matcop(nudim, ncomp, ncomp, nmsh, u, uold) call interp(ncomp, nmsh, xx, nudim, u, * nmold, xxold, uold) endif endif else * err8 is greater than err6 * For a linear problem, set all elements of etest8 to one, * which makes the error test stricter. (The elements of etest8 * may have already been set to one earlier in this routine.) if (linear) call dload(ntol, one, etest8, 1) * Selectively refine the mesh using the old solution and the * 6th order deferred correction. Then, for a nonlinear prpblem, * interpolate the old solution onto the new mesh. ipow = 4 * The array def6 will be overwritten by selmsh. call selmsh(ncomp, nmsh, ntol, ltol, tol, * nfxpnt, fixpnt, ipow, nmax, * xx, ncomp, uold, def6, irefin, ihcomp, * nmold, xxold, ermx, double, maxmsh) if (.not. maxmsh) then if (linear) then call initu(ncomp, nmsh, xx, nudim, u) else call interp(ncomp, nmsh, xx, nudim, u, * nmold, xxold, uold) endif endif * end of logic for err8 greater than err6 endif if (pdebug .and. .not.succes) write(6,903) return 901 format(1h ,'conv8') 902 format(1h ,'err6, err8, er6old, er8old',4(1pe11.3)) 903 format(1h ,'8th order fails error tests.') end subroutine fail8(ncomp, nmsh, nfxpnt, fixpnt, nmax, * ntol, ltol, tol, nmold, * xx, nudim, u, def6, xxold, uold, * ihcomp, irefin, ermx, double, maxmsh) implicit double precision(a-h,o-z) dimension fixpnt(*), ltol(ntol), tol(ntol) dimension xx(*), u(nudim,*), def6(ncomp,*) dimension xxold(*), uold(ncomp,*) dimension ihcomp(*), irefin(*) dimension ermx(*) logical double, maxmsh logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 if (pdebug) write(6,901) * 8th order solution did not converge (the problem must be nonlinear) ipow = 4 * Selectively refine the mesh based on the 6th order deferred * correction and the old solution. * The def6 array is overwritten by selmsh. call selmsh(ncomp, nmsh, ntol, ltol, tol, * nfxpnt, fixpnt, ipow, nmax, * xx, ncomp, uold, def6, irefin, ihcomp, * nmold, xxold, ermx, double, maxmsh) * Interpolate to obtain the new initial solution. if (.not. maxmsh) then call interp(ncomp, nmsh, xx, nudim, u, nmold, xxold, uold) endif return 901 format(1h ,'fail8') end subroutine dccal( ncomp, nmsh, ntol, ltol, * defexp, defimp, dfctol, fval, * ratdc, dfexmx, incmp, inmsh, intol, * derivm, dfimmx, rat1, rat2) implicit double precision (a-h,o-z) dimension ltol(ntol) dimension defexp(ncomp,nmsh-1), defimp(ncomp,nmsh-1) dimension fval(ncomp,nmsh), ratdc(nmsh-1) logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 parameter ( zero = 0.0d+0, one = 1.0d+0 ) parameter ( rtst = 50.0d+0, tstrat = 0.1d+0 ) intrinsic abs, max * blas: idamax * Find dfexmx, the maximum-magnitude element of defexp * in components for which a tolerance is specified. * The component index of dfexmx (incmp), its mesh * interval index (inmsh), and its tolerance index (intol), * are output parameters. dfexmx = zero do 10 it = 1, ntol icmp = ltol(it) idmx = idamax(nmsh-1, defexp(icmp, 1), ncomp) dval = abs(defexp(icmp, idmx)) if (dval .ge. dfexmx) then dfexmx = dval incmp = icmp inmsh = idmx intol = it endif 10 continue if (pdebug) then write(6,901) write(6,902) dfexmx, incmp, inmsh, intol endif * Find derivm (maximum-magnitude element of fval(incmp,*)) * for all mesh points. idmx = idamax(nmsh, fval(incmp, 1), ncomp) derivm = abs(fval(incmp, idmx)) if (pdebug) write(6,903) derivm * For component incmp, go through the mesh intervals to calculate * (1) dfimmx, the maximum implicit deferred correction; * (2) two crucial ratios, rat1 and rat2, used in deciding whether * to refine the mesh; * (3) the array ratdc of deferred-correction ratios (explicit to * implicit). * In defining rat1 and rat2, we consider only intervals for * which the explicit deferred correction (defexp) exceeds the * tolerance dfctol in magnitude. If it does not, the associated * interval does not affect rat1 or rat2, and the value of ratdc * is taken as 1. * rat2 is the maximum-magnitude ratio of sufficiently large * defexp to the larger of defimp and dfctol. * rat1 is the maximum-magnitude ratio of sufficiently large * defexp to the larger in magnitude of defimp and dfctol, but only * for those values of defexp greater than tstrat*dfexmx. * Thus by construction rat1 is less than or equal to rat2. rat1 = zero rat2 = zero dfimmx = zero smtest = tstrat*dfexmx do 100 im = 1, nmsh-1 texp = defexp(incmp, im) timp = defimp(incmp, im) dfimmx = max(dfimmx, abs(timp)) abtexp = abs(texp) if (abtexp .le. dfctol) then ratdc(im) = one else if (abs(timp) .lt. dfctol) timp = dfctol ratdc(im) = texp/timp abrat = abs(ratdc(im)) rat2 = max(rat2, abrat) if (abs(texp) .ge. smtest * .and. abrat .ge. rat1) rat1 = abrat endif * if (pdebug) write(6,904) im, texp, timp, ratdc(im), dfctol 100 continue if (pdebug) write(6,905) rat1, rat2, dfimmx return 901 format(1h ,'dccal') 902 format(1h ,'dfexmx, incmp, inmsh, intol',1pe11.3,3i5) 903 format(1h ,'derivm',1pe11.3) 904 format(1h ,'im, texp, timp, ratdc, dfctol',i5,4(1pe11.3)) 905 format(1h ,'rat1, rat2, dfimmx', 3(1pe11.3)) end subroutine decid4(linear, rat1, rat2, dfexmx, dfimmx, * derivm, dfold, tolval, oldrt1, * onto6, smooth, callrt, strctr, oscchk, double, reposs) implicit double precision (a-h,o-z) logical linear, onto6, smooth, callrt, strctr, * oscchk, double, reposs logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 logical stest parameter ( tenth = 0.1d+0, one = 1.0d+0, two = 2.0d+0 ) parameter ( thrtwo = 32.0d+0 ) parameter ( rtst = 50.0d+0, derval = 50.0d+0 ) * decid4 evaluates information about the deferred corrections * and the nature of the problem, and sets various logical * flags that control the subsequent logic of the algorithm. * This code has been written for clarity, and NOT for efficiency. onto6 = .true. callrt = .false. smooth = .false. oscchk = .false. strctr = .false. reposs = .false. double = .false. * rat2 is always greater than or equal to rat1. if (pdebug) write(6,901) if (pdebug) write(6,902) tolval, rtst stest = .true. if (linear) stest = dfexmx .lt. tenth*dfold if (rat2 .lt. rtst) then if (stest) then smooth = .true. else oscchk = .true. endif return endif * We know now that rat2 .ge. rtst. thttol = thrtwo*tolval if (pdebug) write(6,903) thttol if (rat1 .lt. rtst .and. dfexmx .lt. thttol) then if (stest) then smooth = .true. else oscchk = .true. endif return endif if (rat1 .lt. rtst .and. dfexmx .ge. thttol) then callrt = .true. return endif * We know now that rat1 .ge. rtst (and that rat2 .ge. rtst). if (derivm .gt. derval .and. dfexmx .lt. thttol) then if (stest) then smooth = .true. else oscchk = .true. endif return endif if (derivm .gt. derval .and. dfexmx .gt. thttol) then if (dfimmx .lt. one) then callrt = .true. else strctr = .true. if (linear) then onto6 = .false. if (two*rat1 .ge. oldrt1) double = .true. * end of logic for linear endif * end of logic for dfimmx .ge. one endif return * end of logic for derivm .gt. derval .and dfexmx .gt. thttol endif * To reach this point in the code, both of the following must hold: * rat1 .ge. rtst (which means that rat2 .ge. rtst) * derivm .le. derval * On linear problems, a special termination rule is tried if two * conditions hold: * (1) the 4th order solution has been computed on two consecutive * meshes, one of which is the double of the other * (this holds when double is .true.), and * (2) on both meshes, rat1 .ge. rtst, and derivm is small. When * the conditions in (2) hold for a particular mesh, decid4 * sets reposs to .true. to indicate that Richardson * extrapolation may be possible. * This set of tests is to take care of transients kept out by * initial conditions. if (linear) reposs = .true. return 901 format(1h ,'decid4') 902 format(1h ,'tolval, rtst',2(1pe11.3)) 903 format(1h ,'thttol',1pe11.3) end subroutine dfexcl (ncomp, nmsh, xx, nudim, u, defexp, fval, * tmp, fsub) implicit double precision (a-h, o-z) integer ncomp, nmsh dimension xx(nmsh), u(nudim,nmsh), fval(ncomp, nmsh) dimension defexp(ncomp,nmsh-1), tmp(ncomp,4) external fsub logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 parameter ( half = 0.5d+0, fourth = 0.25d+0, thfrth= 0.75d+0 ) common /consts/ alp1, alp2, alp3, bet0, bet2, bet3, bet4, * a1, b1, c1, d1, e2, f2, c2, d2, b2, a2, * p3, q3, e3, f3, c3, d3, a3, b3, a4, p4, x4, e4, c4, * a5, b5, c5, d5, e5, f5, a6, b6, c6 * Given the nmsh mesh points xx, the estimated solution * u and the array fval of function values at (xx(im), u(*,im)), * im = 1,...,nmsh, dfexcl calculates sixth-order explicit * deferred correction, stored in the array defexp, indexed * over the components and mesh intervals. * The array tmp is workspace for 4 intermediate vectors of * dimension ncomp. do 50 im = 1, nmsh-1 hmsh = xx(im+1) - xx(im) do 10 ic = 1, ncomp tmp(ic,1) = (a5*u(ic, im+1) + b5*u(ic, im)) * + hmsh*(c5*fval(ic,im)- d5*fval(ic,im+1)) tmp(ic,2) = (b5*u(ic,im+1) + a5*u(ic,im)) * + hmsh*(-c5*fval(ic,im+1) + d5*fval(ic,im)) 10 continue call fsub (ncomp, xx(im) + fourth*hmsh, tmp(1,1), * tmp(1,3)) call fsub (ncomp, xx(im) + thfrth*hmsh, tmp(1,2), * tmp(1,4)) do 20 ic = 1, ncomp tmp(ic,1) = half*(u(ic,im+1) + u(ic,im)) * + e5*hmsh*(fval(ic,im+1) - fval(ic,im)) * - f5*hmsh*(tmp(ic,4) - tmp(ic,3)) 20 continue call fsub(ncomp, half*(xx(im) + xx(im+1)), tmp(1,1), * tmp(1,2)) do 30 ic = 1, ncomp defexp(ic,im) = hmsh*(a6*(fval(ic,im+1) + fval(ic,im)) * + b6*(tmp(ic,3) + tmp(ic,4)) + c6*tmp(ic,2)) * - u(ic,im+1) + u(ic,im) 30 continue 50 continue return end subroutine df8cal (ncomp, nmsh, xx, nudim, u, fval, def8, * tmp, fsub) * Given the mesh points xx, the solution u, and the function * values fval, df8cal computes eighth-order deferred corrections, * which are stored in def8. * The array tmp is workspace for 8 intermediate vectors. implicit double precision (a-h, o-z) integer ncomp, nmsh dimension xx(nmsh), u(nudim,nmsh), fval(ncomp,nmsh) dimension def8(ncomp, nmsh-1), tmp(ncomp,8) external fsub logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 parameter ( half = 0.5d+0, two = 2.0d+0 ) parameter ( fc1 = 0.625d+0, fc2= 0.375d+0 ) common /consts/ alp1, alp2, alp3, bet0, bet2, bet3, bet4, * a1, b1, c1, d1, e2, f2, c2, d2, b2, a2, * p3, q3, e3, f3, c3, d3, a3, b3, a4, p4, x4, e4, c4, * a5, b5, c5, d5, e5, f5, a6, b6, c6 do 100 im = 1, nmsh-1 hmsh = xx(im+1) - xx(im) do 10 ic = 1, ncomp tmp(ic,1) = a1*u(ic,im+1) + b1*u(ic,im) * + hmsh*(c1*fval(ic,im+1) + d1*fval(ic,im)) tmp(ic,2) = b1*u(ic,im+1) + a1*u(ic,im) * - hmsh*(c1*fval(ic,im) + d1*fval(ic,im+1)) 10 continue call fsub(ncomp, xx(im) + fc1*hmsh, tmp(1,1), tmp(1,3)) call fsub(ncomp, xx(im) + fc2*hmsh, tmp(1,2), tmp(1,4)) do 20 ic = 1, ncomp tmp(ic,1) = a2*u(ic,im+1) + b2*u(ic,im) * + hmsh*(c2*fval(ic,im+1) + d2*fval(ic,im) * + e2*tmp(ic,3) + f2*tmp(ic,4)) tmp(ic,2) = b2*u(ic,im+1) + a2*u(ic,im) * - hmsh*(d2*fval(ic,im+1) + c2*fval(ic,im) * + f2*tmp(ic,3) + e2*tmp(ic,4)) 20 continue call fsub(ncomp, xx(im) + (half + alp2)*hmsh, tmp(1,1), * tmp(1,5)) call fsub(ncomp, xx(im) + (half - alp2)*hmsh, tmp(1,2), * tmp(1,6)) do 30 ic = 1, ncomp tmp(ic,1) = a3*u(ic,im+1) + b3*u(ic,im) * + hmsh*(c3*fval(ic,im+1) + d3*fval(ic,im) * + e3*tmp(ic,3) + f3*tmp(ic,4) * + p3*tmp(ic,5) + q3*tmp(ic,6)) tmp(ic,2) = b3*u(ic,im+1) + a3*u(ic,im) * - hmsh*(d3*fval(ic,im+1) + c3*fval(ic,im) * + f3*tmp(ic,3) + e3*tmp(ic,4) * + q3*tmp(ic,5) + p3*tmp(ic,6)) 30 continue call fsub (ncomp, xx(im) + (half + alp3)*hmsh, tmp(1,1), * tmp(1,7)) call fsub (ncomp, xx(im) + (half - alp3)*hmsh, tmp(1,2), * tmp(1,8)) do 40 ic = 1, ncomp tmp(ic,1) = a4*(u(ic,im+1) + u(ic,im)) * + hmsh*(c4*(fval(ic,im+1) - fval(ic,im)) * + e4*(tmp(ic,3) - tmp(ic,4)) * + x4*(tmp(ic,7) - tmp(ic,8))) 40 continue call fsub (ncomp, xx(im) + half*hmsh, tmp(1,1), tmp(1,2)) do 50 ic = 1, ncomp def8(ic,im) = * hmsh*(bet0*(fval(ic,im) + fval(ic,im+1)) * + bet2*(tmp(ic,5) + tmp(ic,6)) * + bet3*(tmp(ic,7) + tmp(ic,8)) * + two*bet4*tmp(ic,2)) * - u(ic,im+1) + u(ic,im) 50 continue 100 continue return end subroutine dfimcl( ncomp, nmsh, defexp, chold, dsq, dexr, * ipivot, defimp) implicit double precision (a-h,o-z) dimension defexp(ncomp, nmsh-1), chold(ncomp, ncomp, nmsh-1) dimension dsq(ncomp, ncomp), dexr(ncomp) dimension ipivot(ncomp), defimp(ncomp, nmsh-1) logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 * blas: dcopy parameter (zero = 0.0d+0) * dfimcl calculates the rational deferred correction array, * which is indexed over the components and mesh intervals. call mtload(ncomp, nmsh-1, zero, ncomp, defimp) do 100 im = 1, nmsh-1 call dcopy(ncomp, defexp(1,im), 1, dexr(1), 1 ) do 50 ic = 1, ncomp call dcopy(ncomp, chold(1,ic,im), 1, dsq(1,ic), 1) 50 continue call lufac (ncomp, ncomp, dsq, ipivot, ierlu) if (ierlu .eq. 0) then call lusol (ncomp, ncomp, dsq, ipivot, dexr, * defimp(1, im)) endif 100 continue return end subroutine osc (ncomp, nmsh, dfexmx, incmp, * defcor, ratdc, double, inmsh, onto6, trst6, smooth) implicit double precision (a-h, o-z) dimension defcor(ncomp, *), ratdc(*) logical double, onto6, trst6, smooth logical pdebug common/algprs/ nminit, pdebug, iprint, idum, uval0 intrinsic abs parameter ( zero = 0.0d+0, half = 0.5d+0 ) parameter ( frac1 = 0.1d+0, frac2 = 1.0d-2 ) * For linear problems, subroutine osc performs heuristic tests * to detect an oscillating solution. the tests check whether * (significant) explicit and implicit deferred corrections have * different (componentwise) signs. * dfexmx is the maximum-magnitude explicit deferred correction, * and is known to occur in component incmp. * The array defcor contains the explicit deferred corrections. * The array ratdc contains the ratios of explicit to implicit * deferred corrections. * jsndif counts the number of differences in sign. jsndif = 0 rmax = zero * allsum is the sum of the magnitudes of all deferred corrections, * smlsum is the sum of the magnitudes of the small deferred * corrections, and bigsum is the sum of the magnitudes of the * large deferred corrections. Here, small is defined as less * than half of the maximum. ninter = nmsh - 1 if (pdebug) write(6,901) allsum = zero smlsum = zero bigsum = zero ibig = 0 ism = 0 do 30 im = 1, ninter abdef = abs(defcor(incmp,im)) allsum = allsum + abdef if (abdef .lt. half*dfexmx) then ism = ism + 1 smlsum = smlsum + abdef else ibig = ibig + 1 bigsum = bigsum + abdef endif * The counter of sign differences is incremented if (1) ratdc is negative * (which means that the two deferred corrections have opposite * sign) and (2) the explicit deferred correction is not too small * relative to the maximum. if (pdebug) write(6,902) im, ratdc(im), abdef, frac2*dfexmx if (ratdc(im).lt.zero .and. abdef.ge.frac2*dfexmx) then jsndif = jsndif + 1 * If more than 4 sign differences have occurred, exit after setting * double to .true., which signals that the mesh * should be doubled (i.e., twice as many intervals). if (jsndif.gt.4) then onto6 = .false. double = .true. return endif if (abs(ratdc(im)).ge.rmax) then rmax = abs(ratdc(im)) inmsh = im endif endif 30 continue if (pdebug) write(6,903) rmax, jsndif avsm = zero if (ism.gt.0) avsm = smlsum/ism avbg = zero if (ibig.gt.0) avbg = bigsum/ibig ave = allsum/ninter if (pdebug) write(6,904) ave, avsm, avbg if (avsm.gt.frac1*avbg .or. ave.gt.half*avbg) then * The error appears to be uniformly large. * Signal that the 6th order solution should be calculated. onto6 = .true. elseif (jsndif.eq.0) then * If there were no sign changes, the problem appears to be smooth. smooth = .true. onto6 = .true. else * If the sign changed at between 1 and 4 points, don't go on to * 6th order, and don't ever accept a 6th order solution even if the * error estimate at a later stage indicates that it is OK to do so. * Set double to .false., to signal that the mesh will not necessarily * be doubled. double = .false. onto6 = .false. trst6 = .false. endif return 901 format(1h ,'osc') 902 format(1h ,'im, ratdc, abdef, val',i5,3(1pe11.3)) 903 format(1h ,'rmax, jsndif', 1pe11.3,i5) 904 format(1h ,'ave, avsm, avbg', 3(1pe11.3)) end subroutine ratcor ( ncomp, nmsh, xx, defimp, bhold, dfrat) implicit double precision (a-h, o-z) dimension xx(nmsh), defimp(ncomp,nmsh-1) dimension dfrat(ncomp,nmsh-1), bhold(ncomp,ncomp,nmsh-1) * blas: ddot, dscal parameter (zero = 0.0d+0, half = 0.5d+0, one = 1.0d+0) ninter = nmsh - 1 do 10 im = 1, ninter hmsh = xx(im+1) - xx(im) call dscal(ncomp*ncomp, (-half*hmsh), bhold(1,1,im), 1) 10 continue do 20 im = 1, ninter do 20 ic = 1, ncomp bhold(ic,ic,im) = bhold(ic,ic,im) + one 20 continue do 30 im = 1, ninter do 30 ic = 1, ncomp dfrat(ic,im) = ddot(ncomp, bhold(ic,1,im), ncomp, * defimp(1,im), 1) 30 continue return end subroutine stcons implicit double precision (a-h,o-z) * stcons computes constants needed in integration formulae * and stores them in a labeled common area. parameter ( one = 1.0d+0, four = 4.0d+0, two = 2.0d+0 ) parameter ( five = 5.0d+0, three = 3.0d+0 ) parameter ( half = 0.5d+0, fourth = 0.25d+0