c c Copyright (C) 2002 Vanderbilt group c This file is distributed under the terms of the GNU General Public c License as described in the file 'License' in the current directory. c subroutine pswfnormc(psi,rho,vpot,rcut,l,energy,r, $ rab,dx,con,mesh,ifprt,idim1,idim2) c ----------------------------------- c c recherche des coefficients du polynome implicit none integer l,mesh,ifprt,ndm,idim1,idim2,info real*8 psi(idim1), rho(idim1), vpot(idim1), rcut,energy, + r(idim1), rab(idim1),dx, con(idim2) parameter(ndm=1000) real*8 vpot1(ndm) real*8 r_,rab_,c,c2,amat,y,rc,aenorm integer lam,ipvt,ik common/radmes/r_(ndm),rab_(ndm) common/fundat/ c(6),c2,lam + /funex / amat(6,6),ipvt(6),y(6),rc,aenorm,ik integer stderr,input,iout,ioae,iplot,iologd,iops common /files/ stderr,input,iout,ioae,iplot,iologd,iops real*8 c2l,c2h,precision real*8 funz,rtflsp,chip2,r2 integer i,i1,i2,i3,i4 external funz,rtflsp,chip2 c do i=1,mesh r_(i) = r(i) rab_(i)= rab(i) enddo do i=2,mesh vpot1(i) = vpot(i)/r(i) enddo vpot1(1) = 2.*vpot1(2)-vpot1(3) do i = 1,6 c(i) = 0.0 end do c compute the weight under psi rho(1) = 0.0d0 do 10 i = 2,mesh rho(i) = psi(i)**2 10 continue call radin(mesh,rho,rab,aenorm,idim1) c if (ifprt.gt.-2) write (iout,20) aenorm 20 format (' integral under unpseudised psi is ',f15.10) c do i2 = mesh,1,-2 if ( r(i2) .le. rcut ) goto 110 enddo c 110 continue c c check that the value of i2 is reasonable if ( i2 .lt. 3 .or. i2 .gt. mesh-6 ) then write(iout,*) '***error in subroutine pswfnormc' write(iout,*) 'rcut =',rcut,' generated i2 =',i2, + ' which is illegal with mesh =',mesh call exit(0) endif c rc = r(i2) i1 = i2 - 1 i3 = i2 + 1 i4 = i2 + 2 c rcut = rc if (ifprt.gt.0) then write (iout,*) ' i1, i2, i3, i4, kkbeta =',i1,i2,i3,i4,mesh write (iout,120) ' r1, r2, r3, r4 =',(r(i),i=i1,i4) write (iout,120) 'ps1,ps2,ps3,ps4 =',(psi(i),i=i1,i4) 120 format(a18,4f8.4) write (iout,130) rcut 130 format(' rcut =',f8.4) endif ik = i2 lam = l c aenorm is the norm out to the cutoff aenorm = rho(ik) call fill_matrix(amat,rc,lam) * calcul de mat = lu call dgefa(amat,6,6,ipvt,info) * calcul de y pour ax = y call eval_coeff(r,psi,ik,lam,energy,dx,vpot1,y) c2h = 5.0 ! borne superieure pour c2 c2l = -5.0 ! borne inferieure pour c2 precision = 1e-12 c2 = rtflsp(funz,c2l,c2h,precision) * cherche la valeur de c2 qui verifie l'equ 29a con(1) = c(1) con(2) = c2 do i=2,6 con(i+1)=c(i) enddo c psi(1) = 0.0d0 do i = 2,ik r2 = r(i)*r(i) psi(i) = r(i)**(lam+1) * exp(((((((c(6)*r2+c(5))*r2+c(4))* $ r2+c(3))*r2+c(2))*r2+c2)*r2)+c(1)) enddo do i=1,mesh rho(i) = psi(i)**2 enddo call radin(mesh,rho,rab,aenorm,idim1) if (ifprt.gt.-2) then write (iout,200) aenorm 200 format (' integral under pseudised psi is ',f15.10) endif c test for rhodium return end c c --------------------------------------------------------------- subroutine fill_matrix(amat,rc,l) c -------------------------------- implicit none real*8 amat(6,6),rc integer l c c this routine fills the matrix withn the coefficients taken c from p,p',p'',p''', p(iv), where p is c p(r)= c0 + c4 r^4 + c6 r^6 + c8 r^8 + ... integer pr1(6),cr1(6),pr(6),cr(6),i,j c matrices representing the coefficients (cr) and the powers of r ( pr) data pr1/0,4,6,8,10,12/ data cr1/1,1,1,1,1,1/ do i=1,6 pr(i) = pr1(i) cr(i) = cr1(i) end do do i = 1,5 c fill matrix row by row do j = 1,6 amat(i,j) = cr(j) * rc**(pr(j)*1.0) end do c derivate polynomial expression do j = 1,6 cr(j) = cr(j) * pr(j) pr(j) = max(0,pr(j)-1) end do end do do j = 1,6 amat(6,j) = 0.0 end do amat(6,2) = 2.0*l +5.0 return end c ---------------------------------------------------------------- subroutine eval_coeff (r,psi,ik,l,energy,dx,vpot,y) c ---------------------------------------------------------------- c calcule les coefficients dependant de la fct d'onde calculee c avec tous les electrons. ces coefficients ervent a la resolution c du systeme lineaire. c en entree : ik,nx,vpot comme dans le programme principal c en sortie : une matrice colonne y contenant les coefficients c implicit none integer ik,l,lp1 real*8 r(ik+3),psi(ik+3),vpot(ik+3),energy,dx,y(6) real*8 p,dp,d,vae,dvae,ddvae,rc,rc2,rc3 real*8 deriv_7pts,deriv2_7pts external deriv_7pts,deriv2_7pts,polyn2 real*8 val,d1,d2 c c evaluer p et p' ( dp ) rc = r(ik) c p = psi(ik) c dp = deriv_7pts(psi,ik,rc,dx) call polyn2(rc,val,r(ik-6),psi(ik-6),-1) call polyn2(rc,val,r(ik-6),psi(ik-6), 0) call polyn2(rc, d1,r(ik-6),psi(ik-6),+1) c write(*,*) 'dpsi',p,dp,val,d1 c write(*,*) 'psi',(psi(i),i=ik-1,ik+6) p=val dp = d1 if (p.lt.0.0) then p = -p dp = -dp endif d = dp /p c evaluer vae, dvae, ddvae c vae = vpot(ik) c dvae = deriv_7pts(vpot,ik,rc,dx) c ddvae = deriv2_7pts(vpot,ik,rc,dx) call polyn2(rc,val,r(ik-6),vpot(ik-6),-1) call polyn2(rc,val,r(ik-6),vpot(ik-6), 0) call polyn2(rc, d1,r(ik-6),vpot(ik-6),+1) call polyn2(rc, d2,r(ik-6),vpot(ik-6),+2) write(*,*) 'vpd',vae,val,dvae,d1,ddvae,d2 vae = val dvae= d1 ddvae= d2 c calcul de parametres intervenant dans les calculs successifs lp1 = l + 1 rc2 = rc * rc rc3 = rc2* rc c y(1) = log ( p / rc**lp1 ) y(2) = dp/p - (lp1 / rc) y(3) = vae - energy + (lp1*lp1)/rc2 - d*d y(4) = dvae - 2*(vae - energy + l*lp1/rc2)*d - 2*(lp1*lp1)/rc3 & + 2*(d*d*d) y(5) = ddvae - 2*(dvae - 2*l*lp1/rc3)*d + 6*lp1*lp1/(rc3*rc) & - 2*(vae - energy + l*lp1/rc2 - 3*d*d)* & (vae - energy + l*lp1/rc2 - d*d) return end c ---------------------------------------------------------------- function deriv2_7pts(f,ik,rc,h) c ---------------------------------------------------------------- c c evaluates the second derivative of function f, the function c is given numerically on logarithmic mesh r. c nm = dimension of mesh c ik = integer : position of the point in which the derivative c will be evaluated. c h is distance between x(i) and x(i+1) where c r(i) = exp(x(i))/zed & r(j) = exp(x(j))/zed c implicit none integer a(7),i,ik real*8 f(ik+3),rc,h,sum,sum1,deriv_7pts,deriv2_7pts c integer a(7),i,ik c coefficients for the formula in abramowitz & stegun p.914 c the formula is used for 7 points. data a/4,-54,540,-980,540,-54,4/ ! these are coefficients c formula for linear mesh sum = 0.0 do i=1,7 sum = sum + a(i)*f(i-4+ik) end do sum = 2.0*sum/(720.0*h**2) c transform to logarithmic mesh sum1 = deriv_7pts(f,ik,rc,h) deriv2_7pts = sum/(rc*rc) - sum1 /rc return end c --------------------------------------------------------------- function deriv_7pts(f,ik,rc,h) c --------------------------------------------------------------- c evaluates the first derivative of function f, the function c is given numerically on logarithmic mesh r. c nm = dimension of mesh c ik = integer : position of the point in which the derivative c will be evaluated. c h is distance between x(i) and x(i+1) where c r(i) = exp(x(i))/zed & r(j) = exp(x(j))/zed c implicit none integer a(7),ik,i real*8 f(ik+3),rc,h,sum,deriv_7pts c integer a(7),ik,i c coefficients for the formula in abramowitz & stegun p.914 data a/-12,108,-540,0,540,-108,12/ c formula for linear mesh sum = 0 do i=1,7 sum = sum + a(i)*f(i-4+ik) end do deriv_7pts = sum/(720.0*h) c transform to logarithmic mesh deriv_7pts = deriv_7pts /rc return end c c --------------------------------------------------------------- real*8 function rtflsp(func,x1,x2,xacc) c --------------------------------------------------------------- c using the false position method, find the root of a c function func, known to lie between x1 and x2. c the root, returned as rtflsp, is refined until its accuracy c is +- xacc. c nb.: routine taken from numerical recipes, fortran edition c p 248. c implicit none real*8 func,x1,x2,xacc,fl,fh,zzz,zh,swap,f,del,xl,xh,dx integer maxit,ll,j external func parameter ( maxit = 1000 ) fl = func(x1) fh = func(x2) do while (fl*fh.gt.0) print '(/,'' root not bracketed'')' print '('' actual value of xmin xmax '',2f10.3)',x1,x2 zzz = x1 zh = abs(x2 -x1)/10 print '(/,3x,'' x f(x)'')' do ll = 0,9 zzz = zzz + zh fh = func(zzz) print '(3x,2f10.3)',zzz,fh if (fl*fh .lt. .0d0) then x1 = zzz-zh fl = func(x1) x2 = zzz goto 100 endif end do zzz = (x2-x1) x2 = x2+zzz x1 = x1 - zzz fl = func(x1) fh = func(x2) end do 100 if (fl.lt.0) then xl = x1 xh = x2 else xl = x2 xh = x1 swap = fl fl = fh fh = swap endif dx = xh - xl do j = 1,maxit rtflsp = xl + dx*fl/(fl-fh) f = func(rtflsp) if (f.lt.0) then del = xl - rtflsp xl = rtflsp fl = f else del = xh - rtflsp xh = rtflsp fh = f endif dx = xh - xl if (abs(del).lt.xacc.or.f.eq.0) return end do stop ' rtflsp max iter execeeded' end c -------------------------------------------------------- function funz(x) c -------------------------------------------------------- c cette fonction est la fonction de c2 qu'il faut annuler c pour trouver une valeur de c2 qui verifie la conservation c de la charge de coeur. c cette fonction calcule de plus a chaque fois les ci qui c verifient les equations lineaires avec un c2 donne. c en entree : une valeur de c2 donnee dans x c en sortie, la valeur de la fonction pour cette valeur de x c : c'est la fonction qui correspond a l'equation integrale c implicit none integer ndm,lam,ipvt,ik,i real*8 funz, x, c, c2, amat, y, rc, aenorm, r,rab, $ chip2, psnorm parameter(ndm=1000) common/radmes/ r(ndm),rab(ndm) common /fundat/c(6),c2,lam + /funex/ amat(6,6),ipvt(6),y(6),rc,aenorm,ik external chip2,radin real*8 fint(ndm) c resolution du systeme lineaire pour cette valeur de x (=c2) c(1) = y(1) - x*rc**2 ! ajoute les termes en c2 c(2) = y(2) - 2*x*rc c(3) = y(3) - 2*x c(4) = y(4) ! pas de coeff en c2 c(5) = y(5) c(6) = -x*x call dgesl(amat,6,6,ipvt,c,0) ! resoud le systeme c calcul de la norme de la pseudo-fonction d'onde psnorm = 0.0 fint(1) = 0 do i=2,ik fint(i) = chip2(c,x,lam,r(i)) enddo call radin(ik,fint,rab,psnorm,ndm) funz = log( psnorm / aenorm ) c return end c -------------------------------------------------------- function chip2(c,c2,l,r) c -------------------------------------------------------- implicit none real*8 chip2,c(6),c2,r,r2 integer l r2 = r**2 chip2 = r2**(l+1) * exp(2.0* * (((((((c(6)*r2+c(5))*r2+c(4))*r2+c(3))*r2+c(2))*r2+c2)*r2)+c(1))) return end c c ------------------------------------------------------------------- c ---------------------------------------------------------------- function der3an(l,c,c2,rc) c ---------------------------------------------------------------- c c 3rd derivative of r^(l+1) e^p(r) c implicit none integer l real*8 c(6), c2, rc,pr,dexpr,d2expr,d3expr, der3an external pr,dexpr,d2expr,d3expr der3an= (l-1)*l*(l+1)*rc**(l-2) *exp(pr(c,c2,rc)) + $ 3 * l*(l+1)*rc**(l-1) * dexpr(c,c2,rc) + $ 3 * (l+1)*rc**(l ) *d2expr(c,c2,rc) + $ rc**(l+1) *d3expr(c,c2,rc) return end c ---------------------------------------------------------------- function der4an(l,c,c2,rc) c ---------------------------------------------------------------- c c 4rd derivative of r^(l+1) e^p(r) c implicit none integer l real*8 c(6), c2, rc,pr,dexpr,d2expr,d3expr,d4expr,der4an external pr,dexpr,d2expr,d3expr,d4expr der4an = (l-2)*(l-1)*l*(l+1)*rc**(l-3) * exp(pr(c,c2,rc)) + $ 4 * (l-1)*l*(l+1)*rc**(l-2) * dexpr(c,c2,rc) + $ 6 * l*(l+1)*rc**(l-1) * d2expr(c,c2,rc) + $ 4 * (l+1)*rc**(l ) * d3expr(c,c2,rc) + $ rc**(l+1) * d4expr(c,c2,rc) return end c ---------------------------------------------------------------- function dexpr(c,c2,rc) c ---------------------------------------------------------------- c c 1st derivative of e^p(r) c implicit none real*8 c(6), c2, rc, pr, dpr, dexpr external pr, dpr dexpr= exp(pr(c,c2,rc)) * dpr(c,c2,rc) return end c ---------------------------------------------------------------- function d2expr(c,c2,rc) c ---------------------------------------------------------------- c c 2nd derivative of e^p(r) c implicit none real*8 c(6), c2, rc, pr, dpr, d2pr, d2expr external pr, dpr, d2pr d2expr= exp(pr(c,c2,rc)) * ( dpr(c,c2,rc)**2+d2pr(c,c2,rc) ) return end c ---------------------------------------------------------------- function d3expr(c,c2,rc) c ---------------------------------------------------------------- c c 3nd derivative of e^p(r) c implicit none real*8 c(6), c2, rc, pr, dpr, d2pr, d3pr, d3expr external pr, dpr, d2pr, d3pr d3expr= exp(pr(c,c2,rc)) * ( dpr(c,c2,rc)**3 $ + 3*dpr(c,c2,rc)*d2pr(c,c2,rc) $ + d3pr(c,c2,rc) ) return end c ---------------------------------------------------------------- function d4expr(c,c2,rc) c ---------------------------------------------------------------- c c 4th derivative of e^p(r) c implicit none real*8 c(6), c2, rc, pr, dpr, d2pr, d3pr, d4pr, d4expr external pr, dpr, d2pr, d3pr, d4pr d4expr= exp(pr(c,c2,rc)) * ( dpr(c,c2,rc)**4 $ + 6* dpr(c,c2,rc)**2*d2pr(c,c2,rc) $ + 3*d2pr(c,c2,rc)**2 $ + 4* dpr(c,c2,rc)*d3pr(c,c2,rc) $ + d4pr(c,c2,rc) ) return end c c -------------------------------------------------------- function pr(c,c2,x) c -------------------------------------------------------- c cette fonction evalue le polynome dont les coefficients c sont dans c d'apres la forme suivante c p(x) = c(2)*x^4 + c(3) x^6 + c(4) x^8 + c(5) x^10 c +c(6) x^12 + c2*x^2 + c(1) implicit none real*8 c(6), c2, x, y, pr y = x*x pr = (((((c(6)*y+c(5))*y+c(4))*y+c(3))*y+c(2))*y+c2)*y + + c(1) return end c function dpr(c,c2,x) c -------------------------------------------------------- implicit none real*8 c(6),c2,x,dpr dpr = 2*c2*x + 4*c(2)*x**3 + 6*c(3)*x**5 + 8*c(4)*x**7 & + 10*c(5)*x**9 + 12*c(6)*x**11 return end c function d2pr(c,c2,x) c -------------------------------------------------------- implicit none real*8 c(6),c2,x,d2pr d2pr = 2*c2 + 12*c(2)*x**2 + 30*c(3)*x**4 + 56*c(4)*x**6 & + 90*c(5)*x**8 + 132*c(6)*x**10 return end c function d3pr(c,c2,x) c -------------------------------------------------------- implicit none real*8 c(6),c2,x,d3pr d3pr = 24*c(2)*x + 120*c(3)*x**3 + 336*c(4)*x**5 & + 720*c(5)*x**7 + 1320*c(6)*x**9 return end c function d4pr(c,c2,x) c -------------------------------------------------------- implicit none real*8 c(6),c2,x,d4pr d4pr = 24*c(2) + 360*c(3)*x**2 + 1680*c(4)*x**4 & + 5040*c(5)*x**6 + 11880*c(6)*x**8 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 SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB) INTEGER LDA,N,IPVT(1),JOB DOUBLE PRECISION A(LDA,1),B(1) 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 C INTERNAL VARIABLES 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 * ====================================================================== * NIST Guide to Available Math Software. * Fullsource for module DSCAL from package BLAS1. * Retrieved from NETLIB on Fri Jun 21 17:32:50 1996. * ====================================================================== 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 modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision da,dx(*) integer i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.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