PROGRAM nrelMP ****************************************************************************** * * A routine to determine the parameters in a model potential for an atom with * one valence electron * * (A) Input energies are obtained from the Charlotte Moore Tables in alk.in. * (B) The golden-mean minimization routine from Numerical Recipes is used * to determine the single parameter in the model potential that minimizes * the sum of (Ein(i) - Ecalc(i))^2 * (1) The sum of squares of energies is calculated in ssq.f. This routine * accepts a value of the parameter and a set of initial energies. The * routine calls poten.f to find new energies corresponding to this value * of the parameter. It then evaluates the sum of squares of input - new * energies. * (2) The updated energies are calculated in subroutine poten.f which accepts a * value of the parameter, calculates the potential for this parameter, * and solves the EVP for each input state. * (3) The radial Schrodinger equation EVP is solved in master.f * * ****************************************************************************** implicit doubleprecision(a-h,o-z) character*4 label parameter(NDAT=100) common/energy/wexp(NDAT),wth(NDAT),nda(NDAT),lda(NDAT),ndo common/nucchg/iz common/orblab/label external ssr data tol /1.0d-8/ open(unit=4,file='alk.in',form='formatted',status='old') ndo = 0 *** read in nuclear charge read(4,*) iz DO in = 1,NDAT *** read in remaining orbital data read(4,*,END=100) nda(in),lda(in),wexp(in) call enclab(nda(in),lda(in),*901) write(6,1000) label,wexp(in) 1000 format(4x,a4,f15.6) ndo = ndo+1 END DO 100 continue write(6,1010) ' Enter a,b,c: ' 1010 format(/a,$) read(5,*) a,b,c * The routine golden calls ssr == sum of squares of differences * initially with three trial values (a,b,c) of the potential parameter. * It then refines the paramater until the changes are less than tol. * The final value of the parameter is returned as amin * The final value of ssq is returned as sqer sqer = golden(a,b,c,ssr,tol,amin) write(6,1020) amin,sqer 1020 format(' par =',f10.6,' chi2 =',1pe9.1) DO i = 1,ndo call enclab(nda(i),lda(i),*901) write(6,1030) label,wexp(i),wth(i) 1030 format(4x,a,4x,'wexp =',f10.6,5x,'wth =',f10.6) END DO * * The optimized potential is used to predict the value of the binding * energy for a test case * write(6,*) write(6,*) ' Enter n and l and a guess of E for a test state' nn = ndo+1 read(5,*) nda(nn),lda(nn),wth(nn) ndo = nn call poten(amin) stop 901 stop ' error exit' end FUNCTION GOLDEN(AX,BX,CX,F,TOL,XMIN) IMPLICIT DOUBLEPRECISION(A-H,O-Z) PARAMETER (R=.61803399D0,C=.38196602D0) X0=AX X3=CX IF(ABS(CX-BX).GT.ABS(BX-AX))THEN X1=BX X2=BX+C*(CX-BX) ELSE X2=BX X1=BX-C*(BX-AX) ENDIF write(6,*) ' going in 1st time' F1=F(X1) write(6,*) ' f1 =',F1 F2=F(X2) 1 IF(ABS(X3-X0).GT.TOL*(ABS(X1)+ABS(X2)))THEN IF(F2.LT.F1)THEN X0=X1 X1=X2 X2=R*X1+C*X3 F0=F1 F1=F2 F2=F(X2) ELSE X3=X2 X2=X1 X1=R*X2+C*X0 F3=F2 F2=F1 F1=F(X1) ENDIF GOTO 1 ENDIF IF(F1.LT.F2)THEN GOLDEN=F1 XMIN=X1 ELSE GOLDEN=F2 XMIN=X2 ENDIF RETURN END function ssr(par) implicit doubleprecision(a-h,o-z) parameter(NDAT=100) common/energy/wexp(NDAT),wth(NDAT),nda(NDAT),lda(NDAT),ndo do 50 i = 1,ndo wth(i) = wexp(i) 50 continue call poten(par) ssr = 0d0 do 100 i=1,ndo ssr = ssr + (wth(i)-wexp(i))**2 100 continue return end subroutine poten(par) parameter(NGP=500,NDAT=100) implicit doubleprecision(a-h,o-z) character*4 label common/radial/r(NGP),rp(NGP),rpor(NGP),h & /orblab/label common/energy/wexp(NDAT),wth(NDAT),nda(NDAT),lda(NDAT),ndo common/nucchg/iz dimension p(NGP),q(NGP),z(NGP) data init /0/ if(init.eq.0) then r0 = 0.0005d0 h =.03d0 r(1)=0d0 rp(1)=r0 rpor(1)=0d0 do 100 i=2,NGP rp(i) = r0*dexp((i-1)*h) r(i) = rp(i) - r0 rpor(i) = rp(i)/r(i) 100 continue z0 = iz init = 1 endif do 200 i=1,NGP z(i)= z0 - (z0-1)*(1d0 - exp(-r(i)/par)) 200 continue do 500 itt = 1,ndo n = nda(itt) l = lda(itt) eau = wth(itt) call enclab(n,l,*901) write(6,1150) label,eau 1150 format(/' orbital',a,' E(in) =',f16.11) call master(p,q,z,eau,n,l,m,k,*901) write(6,1200) label,eau,m,k 1200 format(' orbital',a,' E(out) =',f16.11,5x, & ' m =',i4,' k =',i2) wth(itt) = eau 500 continue return 901 stop ' error in potent' end subroutine master(p,q,z,eau,n,l,m,k,*) ************************************************************************ * * find radial eigenfunction and eigenvalue for * Schroedinger equation * p(i) = radial wave function (output) * q(i) = (dp/dr)(i) (output) * z(i) = effective charge; V(r) = -z(r)/r (input) * eau = eigenvalue: input is a guess (input) and * output is precise value (output) * n = principal quantum number (input) * l = angular momentum quantum number (input) * m = practical infinity for orbital (output) * k = number of trys at solution <= ntry (output) * (*) is an error exit ************************************************************************ implicit doubleprecision(a-h,o-z) parameter(NGP=500) common/radial/r(NGP),rp(NGP),rpor(NGP),h dimension p(NGP),q(NGP),z(NGP) dimension u(NGP) data del/1.0d-11/,ntry/30/,alr/800.0d0/ eupper=0.0d0 elower=0.0d0 more=0 less=0 *** nint = number of internal nodes = n-l-1 nint = n - l - 1 **** k = number of trys at solution <= ntry (ntry = 30) k = 0 100 k = k + 1 if(k.gt.ntry) then k = k - 1 return 1 endif **** seek practical infinity: m * * p(r(m)) = exp(- alam r(m) ) * alam = sqrt( -2 (eau + z(m)/r(m)) ) * asymptotic region: alam r(m) > 40 (alr = 40*40/2) * ==> -2 (eau*r(m) + z(m))*r(m) > 1600 * ==> (eau*r(m) + z(m))*r(m) + 800 < 0 * m = NGP + 1 110 continue m = m - 1 if( ( ( eau*r(m) + z(m) )*r(m) + alr ).lt.0.0d0 ) go to 110 **** seek classical turning point: mtp * * classically forbidden region: eau + z(mtp)/r(mtp) < 0 * ==> eau*r(mtp) + z(mtp) < 0 * mtp = m + 1 120 continue mtp = mtp - 1 if( ( eau*r(mtp) + z(mtp) ).lt.0.0d0 ) go to 120 **** integrate inward from practical infinity to classical turning point call insch(p,q,z,eau,l,m,mtp) **** save p(mtp) ptp = p(mtp) qtp = q(mtp) **** integrate outward from origin to classical turning point call outsch(p,q,z,eau,l,mtp) **** match p(mct) to make solution continuous rat = p(mtp)/ptp qtp = rat*qtp mtp1 = mtp + 1 do 130 i = mtp1,m p(i) = rat*p(i) q(i) = rat*q(i) 130 continue **** find number of zeros in p(r) nzero =0 sp = dsign(1.0d0,p(2)) do 140 i = 3,m if(dsign(1.0d0,p(i)).ne.sp) then sp = -sp nzero = nzero + 1 endif 140 continue **** compare nzero with nint * * nzero > nint ==> eau too high, * decrease and try again * nzero < nint ==> eau too low, * increase and try again * nzero = nint ==> eau in correct correct range, * use perturbation theory * if(nzero.gt.nint) then more = more + 1 *** record upper bound on eau if((more.eq.1).or.(eau.lt.eupper)) eupper = eau if(less.ne.0) then **** use average of upper and lower bounds eau = 0.5d0*(eupper+elower) go to 100 endif **** otherwise take 10% less than upper bound eau = 1.10d0*eupper go to 100 elseif(nzero.lt.nint) then less = less + 1 if((less.eq.1).or.(eau.gt.elower)) elower = eau if(more.ne.0) then **** use average of upper and lower bounds eau = 0.5d0*(eupper+elower) go to 100 endif **** otherwise take 10% more than lower bound eau = 0.90d0*elower go to 100 endif **** calculate normalization integral do 150 i = 1,m u(i) = p(i)*p(i)*rp(i) 150 continue anorm = rint(u,1,m,7,h) **** use perturbation theory to calculate energy change * * de = p(mtp) [ q(mpt-) - q(mtp+) ] / [ 2 int( p(r)**2 ) ] * de = p(mtp)*(q(mtp) - qtp)/(2.0d0*anorm) eau = eau + de if((less.ne.0).and.(eau.lt.elower)) then eau = 0.5d0*(eau - de + elower) go to 100 elseif(eau.gt.eupper) then eau = 0.5d0*(eau - de + eupper) go to 100 elseif(dabs(de/eau).gt.del) then go to 100 endif eau = eau - de an = 1.0d0/dsqrt(anorm) do 200 i = 1,m p(i) = an*p(i) q(i) = an*q(i) 200 continue end subroutine adams(p,q,z,eau,l,na,nb) ************************************************************************ * Subroutine to solve the radial Schroedinger equation using * Adams-Moulton NO-point method: NO = 6,7,8,9 * used in conjunction with * subroutines outsch and insch * * p(i) = radial wave function * q(i) = (dp/dr)(i) * dp/dr = q(r) * dq/dr = -2(eau + z(r)/r - l(l+1)/r**2) p(r) * if nb > na then * p(na-NO)...p(na-1) given on input * q(na-NO)...q(na-1) given on input * p(na)...p(nb) found by forward integration * q(na)...q(nb) found by forward integration * else * p(nb+1)...p(nb+NO) given on input * q(nb+1)...q(nb+NO) given on input * p(nb)...p(na) found by backward integration * q(nb)...q(na) found by backward integration * * z(i) effective charge, V(r) = - z(r)/r * eau = energy in a.u. * l = orbital angular momentum quantum number * ************************************************************************ implicit doubleprecision(a-h,o-z) parameter(NGP=500,NO=8) common/radial/r(NGP),rp(NGP),rpor(NGP),h dimension p(NGP),q(NGP),z(NGP) dimension dp(NO),dq(NO),aa(NO) dimension ia(NO) ******** six-point adams method *************************************** * parameter(NO=5) * data ia / 27,-173, 482,-798,1427/ * data id /1440/,iaa/ 475/ ******** seven-point adams method ************************************** * parameter(NO=6) * data ia / -863, 6312,-20211, 37504,-46461, 65112/ * data id /60480/, iaa/19087/ ******** eight-point adams method ************************************** * parameter(NO=7) * data ia / 1375, -11351, 41499, -88547, 123133,-121797, 139849/ * data id / 120960/, iaa/ 36799/ ********** nine-point adams method ************************************* * parameter(NO=8) data ia/ -33953, 312874,-1291214, 3146338,-5033120, 5595358, 1 -4604594, 4467094/ data id /3628800/, iaa /1070017/ ************************************************************************ cof = h/dfloat(id) ang = 0.5d0*l*(l+1) * fill in the preliminary arrays for derivatives if(nb.gt.na) then inc = 1 mstep = nb-na+1 else inc = -1 mstep = na-nb+1 endif i = na-inc*(NO+1) do 100 k = 1,NO i = i + inc dp(k) = inc*rp(i)*q(i) dq(k) = -2*inc*( eau*rp(i) + & (z(i)-ang/r(i))*rpor(i) )*p(i) aa(k) = cof*ia(k) 100 continue a0 = cof*iaa i = na - inc do 500 ii = 1,mstep i = i + inc dpi = inc*rp(i) dqi = -2*inc*( eau*rp(i)+ (z(i)-ang/r(i))*rpor(i) ) b = a0*dpi c = a0*dqi det = 1.0d0 - b*c sp = p(i-inc) sq = q(i-inc) do 200 k = 1,NO sp = sp + aa(k)*dp(k) sq = sq + aa(k)*dq(k) 200 continue p(i) = (sp + b*sq)/det q(i) = (c*sp + sq)/det do 300 k = 1,NO-1 dp(k) = dp(k+1) dq(k) = dq(k+1) 300 continue dp(NO) = dpi*q(i) dq(NO) = dqi*p(i) 500 continue return end subroutine insch(p,q,z,eau,l,nb,na) ************************************************************************ * routine to integrate the radial Schroedinger equation * inward from r(nb) = "practical infinity" to r(na) * * p(i) = radial wave function * q(i) = (dp/dr)(i) * dp/dr = q(r) * dq/dr = -2(eau + z(r)/r - l(l+1)/r**2) p(r) * * z(i) = effective charge, v(r) = -z(r)/r * eau = energy in a.u. * l = orbital angular momentum quantum number * na : r(na) = final point for integration * nb : r(nb) = initial point = practical infinity * nb > na * * Method: * * a) Asymptotic exspansion * p(r) = r**sig exp(-lam r ) [ a0 + a1/r + a2/r**2 + ... * * p(nb)...p(nb-NO) found this way * q(nb)...q(nb-NO) found this way * * b) Continue solution using Adams-Moulton interpolation scheme * * p(nb-NO-1)..p(na) found this way * q(nb-NO-1)..q(na) found this way * ************************************************************************ implicit doubleprecision(a-h,o-z) parameter(NGP=500,NO=8,NX=15) common/radial/r(NGP),rp(NGP),rpor(NGP),h dimension p(NGP),q(NGP),z(NGP) dimension ax(0:NX),bx(0:NX) data eps /1.0d-8/, epr/1d-3/ ************************************************************************ * a) Asymptotic series * * p(r) = r**sig exp(-lam r) [ 1 + a1/r + a2/r**2 + ... ] * lam = sqrt(-2 e) * sig = z/lam * a0 = 1 * a(k+1) = [l(l+1)-(sigk)(sig-k-1)] a(k) /[2 lam (k+1) ] * ************************************************************************ alam = dsqrt(-2.0d0*eau) sig = z(nb)/alam ang = l*(l+1) ax(0) = 1.0d0 bx(0) = -alam do 100 k = 1,NX ax(k) = (ang - (sig-k+1)*(sig-k))*ax(k-1)/(2.0d0*k*alam) bx(k) = ((sig-k+1)*(sig+k) - ang)*ax(k-1)/(2.0d0*k) 100 continue do 300 i = nb,nb-NO,-1 rfac = r(i)**sig*dexp(-alam*r(i)) ps = ax(0) qs = bx(0) rk = 1.0d0 do 200 k = 1,NX rk = rk*r(i) pt = ax(k)/rk qt = bx(k)/rk ps = ps + pt qs = qs + qt xe = dmax1(dabs(pt/ps),dabs(qt/qs)) if(xe.lt.eps) go to 250 200 continue if(xe.gt.epr) write(6,1000) 1000 format(' Warning: asymptotic series not converging ') 250 continue p(i) = rfac*ps q(i) = rfac*qs 300 continue ************************************************************************ * * b) Continue solution inward using Adams method * n.b. be sure NO is consistent in the two routines ! * ************************************************************************ if((nb-NO-1).ge.na) then call adams(p,q,z,eau,l,nb-NO-1,na) endif return end subroutine outsch(p,q,z,eau,l,nb) ************************************************************************ * routine to integrate the radial Schroedinger equation * outward from r(1) = 0 to r(nb) * * p(i) = radial wave function * q(i) = (dp/dr)(i) * dp/dr = q(r) * dq/dr = -2(eau + z(r)/r - l(l+1)/r**2) p(r) * * z(i) = effective charge, v(r) = -z(r)/r * eau = energy in a.u. * l = orbital angular momentum quantum number * nb : r(nb) = final point for integration * * Method: * * a) Factor r**(l+1) and start integration using NO-order Lagrangian * differentiation formulas NO = 6,7,8,9 for LO loops * * p(1)...p(LO*NO+1) found this way * q(1)...q(LO*NO+1) found this way * * b) Continue solution using Adams-Moulton interpolation scheme * * p(LO*NO+2)..p(nb) found this way * q(LO*NO+2)..q(nb) found this way * ************************************************************************ implicit doubleprecision(a-h,o-z) parameter(NGP=500,NO=8,LO=3) common/radial/r(NGP),rp(NGP),rpor(NGP),h dimension p(NGP),q(NGP),z(NGP) dimension us(NO),vs(NO),b(NO),c(NO),d(NO),s(NO) dimension em(NO,NO),fm(NO,NO) dimension ia(NO),ie(NO,NO),ipvt(NO) dimension work(NO),det(2) * * the following coefficients are suitable for a sixth-order method * * use with six-point adams method, set NO = 5 * * data ie / -65, -30, 15, -20, 75, * 1 120, -20, -60, 60, -200, * 2 -60, 60, 20, -120, 300, * 3 20, -15, 30, 65, -300, * 4 -3, 2, -3, 12, 137/ * data ia / -12, 3, -2, 3, -12/ * data id / 60/ * * the following coefficients are suitable for a seventh-order method * * use with seven-point adams method, set NO = 6 * * data ie / -77, -24, 9, -8, 15, -72, * 1 150, -35, -45, 30, -50, 225, * 2 -100, 80, 0, -80, 100, -400, * 3 50, -30, 45, 35, -150, 450, * 4 -15, 8, -9, 24, 77, -360, * 5 2, -1, 1, -2, 10, 147/ / * data ia / -10, 2, -1, 1, -2, 10/ * data id / 60/ * * the following coefficients are suitable for a eighth-order method* * use with eight-point adams method, set NO = 7 * * data ie / -609, -140, 42, -28, 35, -84, 490, * 1 1260, -329, -252, 126, -140, 315, -1764, * 2 -1050, 700, -105, -420, 350, -700, 3675, * 3 700, -350, 420, 105, -700, 1050, -4900, * 4 -315, 140, -126, 252, 329, -1260, 4410, * 5 84, -35, 28, -42, 140, 609, -2940, * 6 -10, 4, -3, 4, -10, 60, 1089/ * data ia / -60, 10, -4, 3, -4, 10, -60/ * data id / 420/ * * the following coefficients are suitable for a ninth-order method* * use with nine-point adams method, set NO = 8 * data ie / -1338, -240, 60, -32, 30, -48, 140, -960, 1 2940, -798, -420, 168, -140, 210, -588, 3920, 2 -2940, 1680, -378, -672, 420, -560, 1470, -9408, 3 2450, -1050, 1050, 0, -1050, 1050, -2450, 14700, 4 -1470, 560, -420, 672, 378, -1680, 2940,-15680, 5 588, -210, 140, -168, 420, 798, -2940, 11760, 6 -140, 48, -30, 32, -60, 240, 1338, -6720, 7 15, -5, 3, -3, 5, -15, 105, 2283/ data ia / -105, 15, -5, 3, -3, 5, -15, 105/ data id / 840/ * data job /01/, nd /NO/, lda /NO/ ************************************************************************ * a) Lagrangian differentiation formulas * * p(r) = r**(l+1) u(r) * * du/dr = v(r) * dv/dr = -2(l+1)/r v(r) -2[ e + z/r ] u(r) * * u = 1 + a1 r + a2 r**2 + ... * v = -z/(l+1) + b1 r + b2 r**2 + ... * ************************************************************************ u0 = 1.0d0 v0 = -z(1)/(l+1) p(1) = 0.0 if(l.eq.0) then q(1) = 1.0d0 else q(1) = 0.0d0 endif do 550 nl = 1,LO i0 = 1 + (nl-1)*NO do 150 i = 1,NO b(i) = id*h*rp(i+i0) c(i) = -2*id*h*(eau*rp(i+i0) + z(i+i0)*rpor(i+i0)) d(i) = -2*id*h*(l+1)*rpor(i+i0) do 100 j = 1,NO if(j.eq.i) then em(i,j) = ie(i,j) - d(i) else em(i,j) = ie(i,j) endif 100 continue 150 continue call dgefa(em,lda,nd,ipvt,info) call dgedi(em,lda,nd,ipvt,det,work,job) do 250 i = 1,NO s(i) = -ia(i)*u0 do 200 j = 1,NO fm(i,j) = ie(i,j) - b(i)*em(i,j)*c(j) s(i) = s(i) - b(i)*em(i,j)*ia(j)*v0 200 continue 250 continue call dgefa(fm,lda,nd,ipvt,info) call dgedi(fm,lda,nd,ipvt,det,work,job) do 350 i = 1,NO us(i) = 0.0d0 do 300 j = 1,NO us(i) = us(i) + fm(i,j)*s(j) 300 continue 350 continue do 450 i = 1,NO vs(i) = 0.0d0 do 400 j = 1,NO vs(i) = vs(i) + em(i,j)*(c(j)*us(j)-ia(j)*v0) 400 continue 450 continue do 500 i = 1,NO p(i+i0) = r(i+i0)**(l+1)*us(i) q(i+i0) = r(i+i0)**l*(r(i+i0)*vs(i)+(l+1)*us(i)) 500 continue u0 = us(NO) v0 = vs(NO) 550 continue *********************************************************************** * * b) Continue solution outward using Adams method * n.b. be sure NO is consistent in the two routines ! * ************************************************************************ na = LO*NO + 2 if(nb.gt.na) then call adams(p,q,z,eau,l,na,nb) endif return end subroutine enclab(n,k,*) c c this program encodes the principal and angular quantum numbers 'n' c and 'l' into an orbital label 'lab' c implicit character*1(l) common/orblab/lab(4) dimension l1(10),l2(10),l3(10),l4(2) data l1/' ','1','2','3','4','5','6','7','8','9'/ data l2/'0','1','2','3','4','5','6','7','8','9'/ data l3/'s','p','d','f','g','h','i','j','k','l'/ data l4/' ','*'/ data nx/10/ if(n.lt.1.or.n.gt.99) then write(6,1000) 1000 format('Error in enclab: principal quantum number out of range') return 1 endif n1=n/10+1 n2=n-n1*10+11 lab(1)=l1(n1) lab(2)=l2(n2) n3=k+1 n4=1 if(n3.lt.0.or.n3.gt.min0(n,nx)) then write(6,1010) 1010 format('Error in enclab: angular quantum number out of range') return 1 endif lab(3)=l3(n3) lab(4)=l4(n4) return end function rint (f,na,nb,nq,h) implicit double precision(a-h,o-z) c c this program calculates the integral of the function f from point na c to point nb using a nq points quadrature ( nq is any integer between c 1 and 14 ). h is the grid size. c written by c. c. j. roothaan c dimension c(105),c1(78),c2(27),d(14),f(nb) equivalence (c1(1),c(1)),(c2(1),c(79)) data c1/ a 1.d0, b 2.d0,1.d0, c 23.d0,28.d0,9.d0, d 25.d0,20.d0,31.d0,8.d0, e 1413.d0,1586.d0,1104.d0,1902.d0,475.d0, f 1456.d0,1333.d0,1746.d0,944.d0,1982.d0,459.d0, g 119585.d0,130936.d0,89437.d0,177984.d0,54851.d0,176648.d0, g 36799.d0, h 122175.d0,111080.d0,156451.d0,46912.d0,220509.d0,29336.d0, h 185153.d0, 35584.d0, i 7200319.d0, 7783754.d0,5095890.d0,12489922.d0,-1020160.d0, i16263486.d0, 261166.d0,11532470.d0,2082753.d0, j 7305728.d0, 6767167.d0, 9516362.d0, 1053138.d0,18554050.d0, j-7084288.d0, 20306238.d0,-1471442.d0,11965622.d0, 2034625.d0, k 952327935.d0, 1021256716.d0, 636547389.d0,1942518504.d0, k-1065220914.d0, 3897945600.d0,-2145575886.d0,3373884696.d0, k -454944189.d0, 1637546484.d0, 262747265.d0, l 963053825.d0, 896771060.d0, 1299041091.d0, -196805736.d0, l 3609224754.d0,-3398609664.d0, 6231334350.d0,-3812282136.d0, l 4207237821.d0, -732728564.d0, 1693103359.d0, 257696640.d0/ data c2 / 5206230892907.d0,5551687979302.d0,3283609164916.d0, m 12465244770050.d0,-13155015007785.d0,39022895874876.d0, m-41078125154304.d0,53315213499588.d0,-32865015189975.d0, m 28323664941310.d0,-5605325192308.d0, 9535909891802.d0, m 1382741929621.d0, n 5252701747968.d0, 4920175305323.d0, 7268021504806.d0, n -3009613761932.d0, 28198302087170.d0,-41474518178601.d0, n 76782233435964.d0,-78837462715392.d0, 81634716670404.d0, n-48598072507095.d0, 34616887868158.d0, -7321658717812.d0, n 9821965479386.d0, 1360737653653.d0/ data d/2.d0,2.d0,24.d0,24.d0,1440.d0,1440.d0,120960.d0, a 120960.d0,7257600.d0,7257600.d0,958003200.d0,958003200.d0, b 5230697472000.d0,5230697472000.d0/ a = 0.0d0 l = na m = nb i = nq*(nq+1)/2 do 100 j = 1,nq a = a + c(i)*( f(l) + f(m) ) l = l + 1 m = m - 1 i = i - 1 100 continue a = a/d(nq) do 200 n = l,m a = a + f(n) 200 continue rint = a*h return end **** contains routines dgedi - dgefa - dgesl subroutine dgedi(a,lda,n,ipvt,det,work,job) integer lda,n,ipvt(1),job double precision a(lda,1),det(2),work(1) c c dgedi computes the determinant and inverse of a matrix c using the factors computed by dgeco or dgefa. c c on entry c c a double precision(lda, n) c the output from dgeco or dgefa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c ipvt integer(n) c the pivot vector from dgeco or dgefa. c c work double precision(n) c work vector. contents destroyed. c c job integer c = 11 both determinant and inverse. c = 01 inverse only. c = 10 determinant only. c c on return c c a inverse of original matrix if requested. c otherwise unchanged. c c det double precision(2) c determinant of original matrix if requested. c otherwise not referenced. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. dabs(det(1)) .lt. 10.0 c or det(1) .eq. 0.0 . c c error condition c c a division by zero will occur if the input factor contains c a zero on the diagonal and the inverse is requested. c it will not occur if the subroutines are called correctly c and if dgeco has set rcond .gt. 0.0 or dgefa has set c info .eq. 0 . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas daxpy,dscal,dswap c fortran dabs,mod c c internal variables c double precision t double precision ten integer i,j,k,kb,kp1,l,nm1 c c c compute determinant c if (job/10 .eq. 0) go to 70 det(1) = 1.0d0 det(2) = 0.0d0 ten = 10.0d0 do 50 i = 1, n if (ipvt(i) .ne. i) det(1) = -det(1) det(1) = a(i,i)*det(1) c ...exit if (det(1) .eq. 0.0d0) go to 60 10 if (dabs(det(1)) .ge. 1.0d0) go to 20 det(1) = ten*det(1) det(2) = det(2) - 1.0d0 go to 10 20 continue 30 if (dabs(det(1)) .lt. ten) go to 40 det(1) = det(1)/ten det(2) = det(2) + 1.0d0 go to 30 40 continue 50 continue 60 continue 70 continue c c compute inverse(u) c if (mod(job,10) .eq. 0) go to 150 do 100 k = 1, n a(k,k) = 1.0d0/a(k,k) t = -a(k,k) call dscal(k-1,t,a(1,k),1) kp1 = k + 1 if (n .lt. kp1) go to 90 do 80 j = kp1, n t = a(k,j) a(k,j) = 0.0d0 call daxpy(k,t,a(1,k),1,a(1,j),1) 80 continue 90 continue 100 continue c c form inverse(u)*inverse(l) c nm1 = n - 1 if (nm1 .lt. 1) go to 140 do 130 kb = 1, nm1 k = n - kb kp1 = k + 1 do 110 i = kp1, n work(i) = a(i,k) a(i,k) = 0.0d0 110 continue do 120 j = kp1, n t = work(j) call daxpy(n,t,a(1,j),1,a(1,k),1) 120 continue l = ipvt(k) if (l .ne. k) call dswap(n,a(1,k),1,a(1,l),1) 130 continue 140 continue 150 continue return end subroutine dgesl(a,lda,n,ipvt,b,job) c c dgesl solves the double precision system c a * x = b or trans(a) * x = b c using the factors computed by dgeco or dgefa. c c on entry c c a double precision(lda, n) c the output from dgeco or dgefa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c ipvt integer(n) c the pivot vector from dgeco or dgefa. c c b double precision(n) c the right hand side vector. c c job integer c = 0 to solve a*x = b , c = nonzero to solve trans(a)*x = b where c trans(a) is the transpose. c c on return c c b the solution vector x . c c error condition c c a division by zero will occur if the input factor contains a c zero on the diagonal. technically this indicates singularity c but it is often caused by improper arguments or improper c setting of lda . it will not occur if the subroutines are c called correctly and if dgeco has set rcond .gt. 0.0 c or dgefa has set info .eq. 0 . c c to compute inverse(a) * c where c is a matrix c with p columns c call dgeco(a,lda,n,ipvt,rcond,z) c if (rcond is too small) go to ... c do 10 j = 1, p c call dgesl(a,lda,n,ipvt,c(1,j),0) c 10 continue c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas daxpy,ddot c integer lda,n,ipvt(1),job double precision a(lda,1),b(1) c double precision ddot,t integer k,kb,l,nm1 c nm1 = n - 1 if (job .ne. 0) go to 50 c c job = 0 , solve a * x = b c first solve l*y = b c if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call daxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue c c now solve u*x = y c do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call daxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue c c job = nonzero, solve trans(a) * x = b c first solve trans(u)*y = b c do 60 k = 1, n t = ddot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 60 continue c c now solve trans(l)*x = y c if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end subroutine dgefa(a,lda,n,ipvt,info) integer lda,n,ipvt(1),info double precision a(lda,1) c c dgefa factors a double precision matrix by gaussian elimination. c c dgefa is usually called by dgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for dgeco) = (1 + 9/n)*(time for dgefa) . c c on entry c c a double precision(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that dgesl or dgedi will divide by zero c if called. use rcond in dgeco for a reliable c indication of singularity. c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas daxpy,dscal,idamax c c internal variables c double precision t integer idamax,j,k,kp1,l,nm1 c c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = idamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (a(l,k) .eq. 0.0d0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue c c compute multipliers c t = -1.0d0/a(k,k) call dscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call daxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (a(n,n) .eq. 0.0d0) info = n return end double precision function dasum(n,dx,incx) c c takes the sum of the absolute values. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dtemp integer i,incx,m,mp1,n,nincx c dasum = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dtemp = dtemp + dabs(dx(i)) 10 continue dasum = dtemp return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,6) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dabs(dx(i)) 30 continue if( n .lt. 6 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,6 dtemp = dtemp + dabs(dx(i)) + dabs(dx(i + 1)) + dabs(dx(i + 2)) * + dabs(dx(i + 3)) + dabs(dx(i + 4)) + dabs(dx(i + 5)) 50 continue 60 dasum = dtemp return end subroutine daxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),da integer i,incx,incy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end double precision function ddot(n,dx,incx,dy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c double precision da,dx(1) integer i,incx,m,mp1,n,nincx c if(n.le.0)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end subroutine dswap (n,dx,incx,dy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal one. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dx(ix) dx(ix) = dy(iy) dy(iy) = dtemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 dtemp = dx(i) dx(i) = dy(i) dy(i) = dtemp dtemp = dx(i + 1) dx(i + 1) = dy(i + 1) dy(i + 1) = dtemp dtemp = dx(i + 2) dx(i + 2) = dy(i + 2) dy(i + 2) = dtemp 50 continue return end integer function idamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c double precision dx(1),dmax integer i,incx,ix,n c idamax = 0 if( n .lt. 1 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end