PROGRAM FCF c Frank-Condon factors c maxb ... maximum excitation c M3 ... number of oscillators c nk0 ... numer of k-points c dimension of the m-vector batch IMPLICIT none integer*4 maxb,nk0,M3,maxe,time,mclock,NB0,NP0 parameter (maxb=100,M3=100,nk0=200,maxe=200,NB0=200,NP0=2000) integer*4 maxn,ii,jj,j,k,kk,nk,n,iq,nv(M3),mv(M3),LEXCL, 1N7,m1,jq,ic,ipiv(M3),i,imaxu,INDX(M3),GEXCL,EEXCL,IWR,nlim, 1state(M3),staten(M3),IEX,it,NEXC,NNEXC,sm(maxe),sn(maxe), 1time0,mvb(NB0*M3),mm,units,NP,ifc real*8 pi,p,ank,ak,akfac,tp,stp,air,aii,ai,BOXF,OH,Cwork(M3), 1xkr(maxb,2*nk0+1),xki(maxb,2*nk0+1),q,xq,fik,box,Determ, 1w(M3,M3),ww(M3,M3),U(M3),UU(M3),e,es(M3),c(M3,M3),fcf1,fcf2, 1d(M3),t(M3,M3),work(3*M3*M3),ti(M3,M3),maxu,d1,d2,ai0,fcf3, 1Det1,Det2,TEMP,Boltzmann,ktcm,EDIF,STRENGTH,t1(M3,M3),EN, 1EE,wcm,wev,wnm,zpeg,zpee,Wi(M3,M3),fcf4,fcf5,fcf6, 1akb(NB0),eeb(nb0),om,S(NP0),XMIN,XMAX,DX,y common/ho/p(maxb,maxb),ank(maxb),e(M3) logical lex,ltest,lspec,ltab,LG character*4 so c Qs excited state normal mode coordinates, qs dimensionless c Ql ground state normal mode coordinates, qs dimensionless c Qs,i=wi,j Ql,j + Ui c qs,i=sqrt(esi/elj)wi,j Ql,j +sqrt(esi) Ui =wwik qli + uui c time=mclock() pi=4.0d0*atan(1.0d0) tp=2.0d0*pi stp=dsqrt(tp) ai=0 do 108 i=1,M3 do 108 j=1,M3 108 w(i,j)=0.0d0 do 109 i=1,M3 nv(i)=1 mv(i)=1 u(i)=0.0d0 w(i,i)=1.0d0 e(i)=1.0d0 109 es(i)=1.0d0 time0=300 inquire(file='fg.inp',exist=lex) if(lex)then inquire(file='fe.inp',exist=lex) if(lex)then write(6,*)'fg.inp and fe.inp found' call ALLIGN(N7,N,WW,UU,e,es,WI) endif endif inquire(file='FCF.OPT',exist=lex) ltest=.false. N7=1 ifc=1 LEXCL=20 nlim=6 IWR=1 NK=6 GEXCL=1 EEXCL=2 M1=10 BOXF=1.0d0 TEMP=298.0d0 OH=1.0d-7 EDIF=7.0d0 STRENGTH=1.0d0 lspec=.true. LG=.false. ltab=.true. XMIN=0.3d0 NP=1000 XMAX=2.0d0 DX=0.006d0 c c UNITS: 1(eV) 2(nm) 3(cm-1) UNITS=1 if(lex)then write(6,*)'FCF.OPT found' open(9,file='FCF.OPT') 92 read(9,919,end=9191,err=9191)so 919 format(a4) write(6,919)so if(so.eq.'TEMP')read(9,*)TEMP if(so.eq.'SPEC')read(9,*)lspec if(so.eq.'TABL')read(9,*)ltab if(so.eq.'UNIT')read(9,*)UNITS if(so.eq.'TIME')read(9,*)TIME0 if(so.eq.'IWR ')read(9,*)IWR if(so.eq.'IFC ')read(9,*)ifc if(so.eq.'GAUS')read(9,*)LG if(so.eq.'TEST')read(9,*)ltest if(so.eq.'EDIF')read(9,*)EDIF if(so.eq.'GEXC')read(9,*)GEXCL if(so.eq.'EEXC')read(9,*)EEXCL if(so.eq.'STRE')read(9,*)STRENGTH if(so.eq.'LEXC')read(9,*)LEXCL if(so.eq.'NLIM')read(9,*)nlim if(so.eq.'NK ')read(9,*)NK if(so.eq.'NP ')read(9,*)NP if(so.eq.'XMIN')read(9,*)XMIN if(so.eq.'XMAX')read(9,*)XMAX if(so.eq.'DX ')read(9,*)DX if(so.eq.'BOXF')read(9,*)BOXF if(so.eq.'M1 ')read(9,*)M1 if(so.eq.'N7 ')read(9,*)N7 if(so.eq.'N ')read(9,*)N if(so.eq.'OH ')read(9,*)OH goto 92 9191 close(9) endif ktcm=0.69507d0*TEMP c maximu exc number/Hermite polynomial degree: if(nk.gt.nk0)call report('too many k-points') if(np.gt.np0)call report('too many spectral points') if(LEXCL.gt.maxe)call report('Too many excitations required') maxn=max(EEXCL,GEXCL) open(22,file='g.txt') NK=3 nv(1)=1 mv(1)=2 87871 maxn=1 do 112 i=N7,N if(mv(i).gt.maxn)maxn=mv(i) 112 if(nv(i).gt.maxn)maxn=nv(i) write(6,6045)maxn-1 6045 format('Max n:',i5) if(maxn.gt.maxb)call report('maxn too high') inquire(file='E.TXT',exist=lex) if(lex)then write(6,*)'E.TXT found' open(9,file='E.TXT') do 110 iq=N7,N 110 read(9,*)e(iq),es(iq) close(9) endif inquire(file='U.TXT',exist=lex) if(lex)then write(6,*)'U.TXT found' open(9,file='U.TXT') read(9,*)(uu(iq),iq=N7,N) close(9) endif inquire(file='W.TXT',exist=lex) if(lex)then write(6,*)'W.TXT found' open(9,file='W.TXT') do 114 iq=N7,N 114 read(9,*)(ww(iq,jq),jq=N7,N) close(9) endif maxu=0.0d0 imaxu=1 do 1131 iq=N7,N if(dabs(uu(iq)).gt.maxu)then maxu=dabs(uu(iq)) imaxu=iq endif 1131 continue c check - Duschinsky matrices: do 113 iq=N7,N u(iq)=uu(iq)/dsqrt(es(iq)) do 113 jq=N7,N c (first index-excited) 113 w(iq,jq)=ww(iq,jq)*dsqrt(e(iq)/es(jq)) c check - dimensionless C matrix: do 201 i=N7,N do 201 j=N7,N t(i,j)=0.0d0 if(i.eq.j)t(i,j)=e(i) do 201 k=N7,N 201 t(i,j)=t(i,j)+w(k,i)*es(k)*w(k,j) call inv(t,N-N7+1,M3,ipiv,work,ti) do 202 i=N7,N do 202 j=N7,N c(i,j)= 2*dsqrt(e(i))*ti(i,j)*dsqrt(e(j)) 202 if(i.eq.j)c(i,j)=c(i,j)-1.0d0 c check - dimensionless D -vector: do 203 i=N7,N d(i)=0.0d0 do 203 j=N7,N do 203 k=N7,N 203 d(i)=d(i)-2*dsqrt(e(i))*ti(i,k)*w(k,j)*es(j)*u(j) c check - I(0,0) factor: do 205 i=N7,N do 205 j=N7,N t(i,j)=0.0d0 205 if(i.eq.j)t(i,j)=e(i)*es(i) call DTRM(t,M3,N-N7+1,Det1,INDX,Cwork) do 206 i=N7,N do 206 j=N7,N t(i,j)=0.0d0 if(i.eq.j)t(i,j)=e(i) do 206 k=N7,N 206 t(i,j)=t(i,j)+w(k,i)*es(k)*w(k,j) do 207 i=N7,N do 207 j=N7,N t1(i,j)=0.0d0 do 207 k=N7,N 207 t1(i,j)=t1(i,j)+w(i,k)*t(k,j) call DTRM(t1,M3,N-N7+1,Det2,INDX,Cwork) d1=0.0d0 d2=0.0d0 do 208 i=N7,N d1=d1-0.5d0*u(i)*es(i)*u(i) do 208 j=N7,N do 208 k=N7,N 208 d2=d2+0.5d0*u(i)*es(i)*w(i,j)*ti(j,k)*w(i,k)*u(i) ai0=dsqrt(dsqrt(Det1**(4*(N-N7+1))))/dsqrt(Det2)*exp(d1+d2) c write(6,6009)imaxu,maxu 6009 format(' Max shift: ',i4,f15.4) box=BOXF*(dsqrt(2.0d0*dble(maxn))*3.0d0*sqrt(2.0d0)+maxu) c iniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniini do 105 k=1,maxn akfac=1.0d0 do 101 kk=1,k-1 101 akfac=akfac*dble(kk) 105 ank(k)=dsqrt(dsqrt(1.0d0/pi)/2.0d0**(k-1)/akfac) do 103 ii=1,maxn do 103 jj=1,maxn 103 p(jj,ii)=0.0d0 p(1,1)=1.0d0 k=1 102 k=k+1 p(k,1 )=-p(k-1,2) do 104 j=2,k-2 104 p(k,j )=(2.0d0*p(k-1,j-1)-dble(j)*p(k-1,j+1)) if(k-2.gt.0)p(k,k-1)=2.0d0*p(k-1,k-2) p(k,k )=2.0d0*p(k-1,k-1) if(k.lt.maxn)goto 102 c iniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniiniini c c transformation of all basis:x(q) .. harmonic oscillator fcion of q c xk=int(x(q) exp(-ik.q),q=-inf...inf) do 1 k=-nk,nk ak=tp*dble(k)/box do 1 j=1,maxn xkr(j,k+nk+1)=0.0d0 xki(j,k+nk+1)=0.0d0 do 1 ii=1,j c int(q^n exp(-q^2/2-ik.q)/sqrt(2*pi),q=-inf...inf) c air,aii .. real an imag. parts call intq(ii-1,ak,air,aii) xkr(j,k+nk+1)=xkr(j,k+nk+1)+p(j,ii)*ank(j)*air*stp/box 1 xki(j,k+nk+1)=xki(j,k+nk+1)+p(j,ii)*ank(j)*aii*stp/box write(6,*)'trafo done' open(30,file='0.prn') open(39,file='n.prn') open(40,file='0f.prn') open(49,file='nf.prn') do 2 iq=-500,500 q=dble(iq)/1000.d0*box xq=0.0d0 do 3 kk=-nk,nk ak=tp*dble(kk)/box 3 xq=xq+xkr(1,kk+nk+1)*cos(ak*q)-xki(1,kk+nk+1)*sin(ak*q) write(30,3001)q,xq 3001 format(2f20.3) write(40,3001)q,fik(1,q) xq=0.0d0 do 39 kk=-nk,nk ak=tp*dble(kk)/box 39 xq=xq+xkr(maxn,kk+nk+1)*cos(ak*(q-maxu)) 1 -xki(maxn,kk+nk+1)*sin(ak*(q-maxu)) write(39,3001)q,xq 2 write(49,3001)q,fik(maxn,q-maxu) close(30) close(40) close(39) close(49) call DTRM(ww,M3,N-N7+1,Determ,INDX,Cwork) write(6,6111)Determ,N7,N 6111 format(' Transformation determinant:',f10.3,' N7:',i2,' N:',i3) Determ=dsqrt(Determ) zpeg=0.0d0 zpee=0.0d0 do 52 iq=N7,N zpeg=zpeg+e(iq)/2.0d0 52 zpee=zpee+es(iq)/2.0d0 if(ltest)goto 8787 it=0 if(ltab)then open(9,file='FCF.TAB') write(9,*) write(9,*) write(9,*) endif if(lspec)then do 400 iq=1,NP 400 S(iq)=0.0d0 endif c produce ground states c loop start 11111111111111111111111111111111111111111111111111 do 31 Nnexc=0,GEXCL do 411 iex=1,Nnexc 411 sn(iex)=N7 51 DO 19661 I=N7,N 19661 STATEn(I)=0 do 44 iex=1,Nnexc 44 staten(sn(iex))=staten(sn(iex))+1 EN=0.0d0 DO 45 I=N7,N 45 EN=EN+dble(staten(I))*E(I) Boltzmann=exp(-EN/ktcm) if(Boltzmann.lt.OH)goto 444 if(iwr.ne.0)write(6,6112)Boltzmann 6112 format(' exp(-E/kt):',f10.6) c do 45 iex=1,Nnexc c5 if(staten(sn(iex)).gt.m1)goto 444 mm=0 c produce excited states c 22222222222222222222222222222222222222222222222222222222222222 do 311 Nexc=0,EEXCL do 4111 iex=1,Nexc 4111 sm(iex)=N7 511 DO 661 I=N7,N 661 STATE(I)=0 do 441 iex=1,Nexc 441 state(sm(iex))=state(sm(iex))+1 ee=0.0d0 DO 46 I=N7,N 46 ee=ee+dble(state(I))*ES(I) c do 45 iex=1,Nnexc c5 if(staten(sn(iex)).gt.m1)goto 344 do 47 iq=N7,N nv(iq)=staten(iq)+1 47 mv(iq)=state(iq) +1 it=it+1 8787 if(iwr.ne.0)then write(6,6001) 6001 format(/,'<',$) write(6,6002)(nv(iq)-1,iq=N7,N) 6002 format(i1,$) write(6,6005) 6005 format('|',$) write(6,6002)(mv(iq)-1,iq=N7,N) write(6,6004) 6004 format('>',/) endif wcm=8065.5d0*EDIF+ee-en+zpee-zpeg wnm=1.0d7/wcm wev=wcm/8065.5d0 if(units.eq.1)then om=wev else if(units.eq.2)then om=wnm else if(units.eq.3)then om=wcm else call report('unknown units') endif endif endif if(ifc.eq.1)then ak=fcf1(nv,mv,nk,ww,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH) else if(ifc.eq.2)then ak=fcf2(nv,mv,nk,ww,uu,LEXCL,nlim,xkr,xki,N7,N,box, 1 m1,ai,OH) else if(ifc.eq.3)then ak=fcf3(nv,mv,nk,wi,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH) else if(ifc.eq.4)then ak=fcf4(nv,mv,nk,ww,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH) else if(ifc.eq.5)then ak=fcf5(nv,mv,nk,wi,uu,LEXCL,xkr,xki,N7,N,box,m1,ai,OH) else if(ifc.eq.6)then ak=fcf6(nv,mv,nk,wi,uu,xkr,xki,N7,N,box,ai,OH,time0) else if(ifc.eq.7.or.ifc.eq.8)then if(GEXCL.ne.0)call report('ifc 7 and 8 only for ground 0') mm=mm+1 eeb(mm)=om do 301 iq=N7,N 301 mvb(mm+(iq-1)*NB0)=mv(iq) if(mm.eq.NB0)then if(ifc.eq.7)then call fcf7(nv,mvb,nk,wi,uu,LEXCL,xkr,xki,N7,N,box, 1 m1,akb,OH,mm) else call fcf8(nv,mvb,nk,wi,uu,xkr,xki,N7,N,box, 1 akb,OH,time0,mm) endif do 303 ii=1,mm ak=akb(ii) ak=ak*Determ om=eeb(ii) y=STRENGTH*ak**2*Boltzmann if(lspec)call tabprnfs(S,XMIN,XMAX,LG,NP,DX,om,y,0) 303 if(ltab)write(9,900)it,om,y mm=0 endif else call report('unknown IFC') endif endif endif endif endif endif endif if(ifc.ne.7.and.ifc.ne.8)then ak=ak*Determ ai=ai*Determ if(iwr.ne.0)write(6,6003)ifc,ak,ai 6003 format(' FCF',i1,' = ',2f10.6) if(mv(1).lt.8.and.nv(1).eq.1)then mv(1)=mv(1)+1 write(22,222)ak 222 format(f12.6,$) goto 87871 endif if(nv(1).eq.1)then nv(1)=nv(1)+1 mv(1)=3 goto 87871 endif write(22,222)ak if(nk.lt.20)then write(22,221)nk 221 format(i3) nk=nk+1 nv(1)=1 mv(1)=2 goto 87871 endif close(22) if(ltest)stop y=STRENGTH*ak**2*Boltzmann if(lspec)call tabprnfs(S,XMIN,XMAX,LG,NP,DX,om,y,0) if(ltab)write(9,900)it,om,y 900 format(i4,f18.9,g14.6) endif 344 if(Nexc.eq.0)goto 311 do 811 ic=Nexc,1,-1 811 if(sm(ic).lt.N)goto 911 goto 311 911 do 1011 i=ic+1,Nexc 1011 sm(i)=sm(ic)+1 sm(ic)=sm(ic)+1 goto 511 311 continue c loop end 11111111111111111111111111111111111111111111111111 c 22222222222222222222222222222222222222222222222222222222222222 444 if(Nnexc.eq.0)goto 31 do 81 ic=Nnexc,1,-1 81 if(sn(ic).lt.N)goto 91 goto 31 91 do 401 i=ic+1,Nnexc 401 sn(i)=sn(ic)+1 sn(ic)=sn(ic)+1 goto 51 31 continue c loop end 11111111111111111111111111111111111111111111111111 if(ifc.eq.7.or.ifc.eq.8)then if(ifc.eq.7)then call fcf7(nv,mvb,nk,w,u,LEXCL,xkr,xki,N7,N,box, 1 m1,akb,OH,mm) else call fcf8(nv,mvb,nk,w,u,xkr,xki,N7,N,box, 1 akb,OH,time0,mm) endif do 304 ii=1,mm ak=akb(ii) ak=ak*Determ om=eeb(ii) y=STRENGTH*ak**2*Boltzmann if(lspec)call tabprnfs(S,XMIN,XMAX,LG,NP,DX,om,y,0) 304 if(ltab)write(9,900)it,om,y endif if(ltab)write(9,90011) if(ltab)close(9) 90011 format(60(1h-)) if(lspec)call tabprnfs(S,XMIN,XMAX,LG,NP,DX,0.0d0,0.0d0,1) write(6,*)' Time ',mclock()-time stop end c ======================================================= subroutine tabprnfs(S,wsmin,wsmax,LG,npoint,hw,wo,x,ic) implicit none integer*4 ic,npoint,i real*8 wsmax,wsmin,s(*),x,dw,dd,hw,xx,wo,sf logical LG if(ic.eq.0)then c disperse value over spectrum dw=(wsmax-wsmin)/dble(npoint-1) if(LG)then dd=hw/2.0d0/dsqrt(log(2.0d0)) else dd=hw/2.0d0 endif do 4 i=1,npoint xx=wsmin+dw*dble(i-1) 4 s(i)=s(i)+sf(dd,LG,xx,wo)*x else c write spectrum call wrspec2('S.PRN',npoint,wsmin,dw,s) endif return end c ================================================================== subroutine wrspec2(sn,nw,wmin,dw,spec) implicit none character*(*) sn integer*4 nw,iw real*8 wmin,dw,spec(*),w open(8,file=sn) do 1 iw=1,nw w=wmin+dble(iw-1)*dw 1 write(8,805)w,spec(iw) 805 format(f18.7,g15.6) close(8) write(6,*)sn,' written' return end c ================================================================== function sf(d,g,x,x0) implicit none logical g real*8 pi,spi,sf,d,x,x0,dd pi=4.0d0*atan(1.0d0) spi=dsqrt(pi) dd=((x-x0)/d)**2 c write(6,*)d,x,x0,g if(g)then if(dd.lt.20.0d0)then sf=exp(-dd)/d/spi else sf=0.0d0 endif else if(dd.lt.1.0d20)then sf=1.0d0/d/(dd+1.0d0)/pi else sf=0.0d0 endif endif return end c ------------------------------------------------------------ function fcf1(nv,mv,nk,w,u,LEXCL,xkr,xki,N7,N,box, 1m1,fcsumi,OH) c via delta-f, standard implicit none integer*4 maxb,nk0,M3,maxe parameter (maxb=100,M3=100,nk0=200,maxe=200) integer*4 nv(*),mv(*),N7,i,nk,N,Nexc,iex,m1,j, 1state(M3),si(maxe),kstate(M3),ic,LEXCL,ii real*8 w(M3,M3),u(*),fcf1,fcsumr,fcsumi,xkr(maxb,2*nk0+1), 1xki(maxb,2*nk0+1),prs,pis,prl,pil,box,fim,pi,tp,aw(M3), 1ptr,pti,OH,t,xkrn,xkin,ak,air,aii,ank,p,stp,e common/ho/p(maxb,maxb),ank(maxb),e(M3) c fcf Frank-Condon factor c mv,nv excited and ground quantum numbers (1=ground state,nosc=0) c nk number of k-expansion points (k=-nk ... nk) c qs=wql+u c LEXCL maximum sum of k-numbers c xkr,xki real and imaginary parts of the exp. coefficients c N7..N number of oscillators, typically N7=1,N=3*nat-6 c box=L the NM interval (-L/2...L/2) c m1 maximum excitation on one center pi=4.0d0*atan(1.0d0) tp=2.0d0*pi stp=dsqrt(tp) fcsumr=0.0d0 fcsumi=0.0d0 c distribute k-"excitations" on electronic excited coeff C "Ground state" zero coeffs. - always: c a.b = (ReaReb-ImaImb) + i (ReaImb+ImaReb): c ps=xs1(0)*xs2(0)...xsN(0) prs= xkr(mv(N7),nk+1) pis= xki(mv(N7),nk+1) do 1 i=N7+1,N t = prs*xkr(mv(i),nk+1)-pis*xki(mv(i),nk+1) pis= prs*xki(mv(i),nk+1)+pis*xkr(mv(i),nk+1) prs=t 1 if(abs(prs)+abs(pis).lt.OH)goto 1001 c pl=xl1(0)*xl2(0)...xlN(0) prl= xkr(nv(N7),nk+1) pil= xki(nv(N7),nk+1) do 2 i=N7+1,N t = prl*xkr(nv(i),nk+1)-pil*xki(nv(i),nk+1) pil= prl*xki(nv(i),nk+1)+pil*xkr(nv(i),nk+1) prl=t 2 if(abs(prl)+abs(pil).lt.OH)goto 1001 c ps* x pl: c a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb): fcsumr=fcsumr+(prs*prl+pis*pil)*box fcsumi=fcsumi+(prs*pil-pis*prl)*box c c distribute LEXCL onto the excited k-expansions N7...N c loop start 22222222222222222222222222222222222222222222222222 1001 do 3 Nexc=1,LEXCL c initial indices - all to the first center: do 41 iex=1,Nexc 41 si(iex)=N7 5 DO 1966 I=N7,N 1966 STATE(I)=0 do 4 iex=1,Nexc 4 state(si(iex))=state(si(iex))+1 do 42 iex=1,Nexc 42 if(state(si(iex)).gt.m1)goto 333 c c odd excitation - positive (N+1)/2 c even excitation - neagtive N/2 do 43 I=N7,N if(mod(state(I),2).eq.0)then kstate(I)=-state(I)/2 else kstate(I)=(state(I)+1)/2 endif if(abs(kstate(i)).gt.nk)goto 333 43 continue c c calculate the phase and product exp(if)xs... fim=0.0d0 do 10 I=N7,N 10 fim=fim+u(I)*kstate(I)*tp/box prs=cos(fim) pis=sin(fim) do 11 i=N7,N t = prs*xkr(mv(i),nk+1+kstate(i))-pis*xki(mv(i),nk+1+kstate(i)) pis= prs*xki(mv(i),nk+1+kstate(i))+pis*xkr(mv(i),nk+1+kstate(i)) prs=t 11 if(abs(prs)+abs(pis).lt.OH)goto 333 c write(6,*)'ks',kstate(1),prs,pis c c calculate the aw products (skip 2 pi /L multiplier) do 12 i=N7,N aw(i)=0.0d0 do 12 j=N7,N c w-first index excited 12 aw(i)=aw(i)+dble(kstate(j))*w(j,i)*tp/box c c calculate the phase and product (exp(-if)xs*(1)xs*(2)...)xl(1)... ptr=1.0d0 pti=0.0d0 do 111 i=N7,N ak=aw(i) xkrn=0.0d0 xkin=0.0d0 do 6 ii=1,nv(i) call intq(ii-1,ak,air,aii) xkrn=xkrn+p(nv(i),ii)*ank(nv(i))*air*stp 6 xkin=xkin+p(nv(i),ii)*ank(nv(i))*aii*stp c write(6,*)xkrn,xkin c write(6,*)xkr(nv(i),nk+1+kstate(i)),xki(nv(i),nk+1+kstate(i)) t = ptr*xkrn-pti*xkin pti= ptr*xkin+pti*xkrn ptr=t 111 if(abs(ptr)+abs(pti).lt.OH)goto 333 c write(6,*)'kl',kstate(1),ptr,pti c a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb): fcsumr=fcsumr+(prs*ptr+pis*pti) fcsumi=fcsumi+(prs*pti-pis*ptr) c find index to be changed 333 do 8 ic=Nexc,1,-1 8 if(si(ic).lt.N)goto 9 goto 3 9 do 102 i=ic+1,Nexc 102 si(i)=si(ic)+1 si(ic)=si(ic)+1 goto 5 3 continue c loop end 22222222222222222222222222222222222222222222222222 fcf1=fcsumr return end c ------------------------------------------------------------ function fcf2(nv,mv,nk,w,u,LEXCL,nlim,xkr,xki,N7,N,box, 1m1,fcsumi,OH) c left sin ... very slow, results virtually sae as fcf1 c not recommended for normal computation implicit none integer*4 maxb,nk0,M3,maxe parameter (maxb=100,M3=100,nk0=200,maxe=200) integer*4 nv(*),mv(*),nlim,N7,i,nk,N,Nexc,Nnexc,iex,m1,j, 1state(M3),staten(M3),si(maxe),sn(maxe),kstate(M3),kstaten(M3), 1naw(M3),ic,LEXCL real*8 w(M3,M3),u(*),fcf2,fcsumr,fcsumi,xkr(maxb,2*nk0+1), 1xki(maxb,2*nk0+1),prs,pis,prl,pil,box,fim,pi,tp,aw(M3), 1dk,sf,ptr,pti,OH,t c fcf Frank-Condon factor c mv,nv excited and ground quantum numbers (1=ground state,nosc=0) c nk number of k-expansion points (k=-nk ... nk) c qs=wql+u c LEXCL maximum sum of k-numbers c nlim .. take ground ks only within this limit around excited ks c xkr,xki real and imaginary parts of the exp. coefficients c N7..N number of oscillators, typically N7=1,N=3*nat-6 c box=L the NM interval (-L/2...L/2) c m1 maximum excitation on one center pi=4.0d0*atan(1.0d0) tp=2.0d0*pi fcsumr=0.0d0 fcsumi=0.0d0 c distribute k-"excitations" on electronic excited coeff C "Ground state" zero coeffs. - always: c a.b = (ReaReb-ImaImb) + i (ReaImb+ImaReb): c ps=xs1(0)*xs2(0)...xsN(0) prs= xkr(mv(N7),nk+1) pis= xki(mv(N7),nk+1) do 1 i=N7+1,N t = prs*xkr(mv(i),nk+1)-pis*xki(mv(i),nk+1) pis= prs*xki(mv(i),nk+1)+pis*xkr(mv(i),nk+1) prs=t 1 if(abs(prs)+abs(pis).lt.OH)goto 1001 c pl=xl1(0)*xl2(0)...xlN(0) prl= xkr(nv(N7),nk+1) pil= xki(nv(N7),nk+1) do 2 i=N7+1,N t = prl*xkr(nv(i),nk+1)-pil*xki(nv(i),nk+1) pil= prl*xki(nv(i),nk+1)+pil*xkr(nv(i),nk+1) prl=t 2 if(abs(prl)+abs(pil).lt.OH)goto 1001 c ps* x pl: c a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb): fcsumr=fcsumr+(prs*prl+pis*pil)*box fcsumi=fcsumi+(prs*pil-pis*prl)*box c c distribute LEXCL onto the excited k-expansions N7...N c loop start 22222222222222222222222222222222222222222222222222 1001 do 3 Nexc=1,LEXCL c initial indices - all to the first center: do 41 iex=1,Nexc 41 si(iex)=N7 5 DO 1966 I=N7,N 1966 STATE(I)=0 do 4 iex=1,Nexc 4 state(si(iex))=state(si(iex))+1 do 42 iex=1,Nexc 42 if(state(si(iex)).gt.m1)goto 333 c c odd excitation - positive (N+1)/2 c even excitation - neagtive N/2 do 43 I=N7,N if(mod(state(I),2).eq.0)then kstate(I)=-state(I)/2 else kstate(I)=(state(I)+1)/2 endif if(abs(kstate(i)).gt.nk)goto 333 43 continue c c calculate the phase and product exp(if)xs... fim=0.0d0 do 10 I=N7,N 10 fim=fim+u(I)*kstate(I)*tp/box prs=cos(fim) pis=sin(fim) do 11 i=N7,N t = prs*xkr(mv(i),nk+1+kstate(i))-pis*xki(mv(i),nk+1+kstate(i)) pis= prs*xki(mv(i),nk+1+kstate(i))+pis*xkr(mv(i),nk+1+kstate(i)) prs=t 11 if(abs(prs)+abs(pis).lt.OH)goto 333 c write(6,*)'ks',kstate(1),prs,pis c c calculate the aw products (skip 2 pi /L multiplier) do 12 i=N7,N aw(i)=0.0d0 do 12 j=N7,N c w-first index excited 12 aw(i)=aw(i)+dble(kstate(j))*w(j,i) c closest integer value: do 121 i=N7,N 121 naw(i)=NINT(aw(i)) c c take n-k vectors in the vicinity of aw c generate excitations on the n(ground state) site c loop start 11111111111111111111111111111111111111111111111111 do 31 Nnexc=1,LEXCL do 411 iex=1,Nnexc 411 sn(iex)=N7 51 DO 19661 I=N7,N 19661 STATEn(I)=0 do 44 iex=1,Nnexc 44 staten(sn(iex))=staten(sn(iex))+1 do 45 iex=1,Nnexc 45 if(staten(sn(iex)).gt.m1)goto 444 do 421 I=N7,N if(mod(staten(I),2).eq.0)then kstaten(I)=-staten(I)/2 else kstaten(I)=(staten(I)+1)/2 endif c take only states that are close to aw if(abs(kstaten(i)).gt.nk)goto 444 421 if(iabs(kstaten(I)-naw(I)).gt.nlim)goto 444 c c calculate the phase and product (exp(-if)xs*(1)xs*(2)...)xl(1)... ptr=1.0d0 pti=0.0d0 do 111 i=N7,N t = ptr*xkr(nv(i),nk+1+kstaten(i))-pti*xki(nv(i),nk+1+kstaten(i)) pti= ptr*xki(nv(i),nk+1+kstaten(i))+pti*xkr(nv(i),nk+1+kstaten(i)) ptr=t if(abs(ptr)+abs(pti).lt.OH)goto 444 dk=dble(kstaten(i))-aw(i) if(dabs(dk).gt.OH)then sf=2.0d0*sin(pi*dk)/(tp*dk/box) else sf=box endif ptr=ptr*sf 111 pti=pti*sf c write(6,*)'kl',kstaten(1),ptr,pti c a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb): fcsumr=fcsumr+(prs*ptr+pis*pti) fcsumi=fcsumi+(prs*pti-pis*ptr) c 444 do 81 ic=Nnexc,1,-1 81 if(sn(ic).lt.N)goto 91 goto 31 91 do 101 i=ic+1,Nnexc 101 sn(i)=sn(ic)+1 sn(ic)=sn(ic)+1 goto 51 31 continue c loop end 11111111111111111111111111111111111111111111111111 c find index to be changed 333 do 8 ic=Nexc,1,-1 8 if(si(ic).lt.N)goto 9 goto 3 9 do 102 i=ic+1,Nexc 102 si(i)=si(ic)+1 si(ic)=si(ic)+1 goto 5 3 continue c loop end 22222222222222222222222222222222222222222222222222 fcf2=fcsumr return end c ------------------------------------------------------------ c ------------------------------------------------------------ function fcf3(nv,mv,nk,w,u,LEXCL,xkr,xki,N7,N,box, 1m1,fcsumi,OH) c ground state first, otherwise same as fcf1 implicit none integer*4 maxb,nk0,M3,maxe parameter (maxb=100,M3=100,nk0=200,maxe=200) integer*4 nv(*),mv(*),N7,i,nk,N,Nexc,iex,m1,j, 1state(M3),si(maxe),kstate(M3),ic,LEXCL,ii real*8 w(M3,M3),u(*),fcf3,fcsumr,fcsumi,xkr(maxb,2*nk0+1), 1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp, 1ptr,pti,OH,t,xkrn,xkin,ak,air,aii,ank,p,stp,e common/ho/p(maxb,maxb),ank(maxb),e(M3) c fcf Frank-Condon factor c mv,nv excited and ground quantum numbers (1=ground state,nosc=0) c nk number of k-expansion points (k=-nk ... nk) c qs=wql+u c LEXCL maximum sum of k-numbers c xkr,xki real and imaginary parts of the exp. coefficients c N7..N number of oscillators, typically N7=1,N=3*nat-6 c box=L the NM interval (-L/2...L/2) c m1 maximum excitation on one center pi=4.0d0*atan(1.0d0) tp=2.0d0*pi stp=dsqrt(tp) fcsumr=0.0d0 fcsumi=0.0d0 c distribute k-"excitations" on electronic excited coeff c start with ground as it may need less k-points do 3 Nexc=0,LEXCL c initial indices - all to the first center: do 41 iex=1,Nexc 41 si(iex)=N7 5 DO 1966 I=N7,N 1966 STATE(I)=0 do 4 iex=1,Nexc 4 state(si(iex))=state(si(iex))+1 do 42 iex=1,Nexc 42 if(state(si(iex)).gt.m1)goto 333 c c odd excitation - positive (N+1)/2 c even excitation - neagtive N/2 do 43 I=N7,N if(mod(state(I),2).eq.0)then kstate(I)=-state(I)/2 else kstate(I)=(state(I)+1)/2 endif if(abs(kstate(i)).gt.nk)goto 333 43 continue c c calculate the product xl... prs=1.0d0 pis=0.0d0 do 11 i=N7,N t = prs*xkr(nv(i),nk+1+kstate(i))-pis*xki(nv(i),nk+1+kstate(i)) pis= prs*xki(nv(i),nk+1+kstate(i))+pis*xkr(nv(i),nk+1+kstate(i)) prs=t 11 if(abs(prs)+abs(pis).lt.OH)goto 333 c c calculate the product xs... ptr=1.0d0 pti=0.0d0 fim=0.0d0 do 111 i=N7,N c calculate corresponding ak: ak=0.0d0 do 12 j=N7,N c w-first index excited 12 ak=ak+dble(kstate(j))*w(i,j) ak=ak*tp/box c calculate the phase: fim=fim+u(i)*ak c c calculate Fourier expansion of excited at ak: xkrn=0.0d0 xkin=0.0d0 do 6 ii=1,mv(i) call intq(ii-1,ak,air,aii) xkrn=xkrn+p(mv(i),ii)*ank(mv(i))*air*stp 6 xkin=xkin+p(mv(i),ii)*ank(mv(i))*aii*stp t = ptr*xkrn-pti*xkin pti= ptr*xkin+pti*xkrn ptr=t 111 if(abs(ptr)+abs(pti).lt.OH)goto 333 t = ptr*cos(fim)-pti*sin(fim) pti= ptr*sin(fim)+pti*cos(fim) ptr=t c a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb): fcsumr=fcsumr+(prs*ptr+pis*pti) fcsumi=fcsumi+(prs*pti-pis*ptr) c find index to be changed 333 if(Nexc.eq.0)goto 3 do 8 ic=Nexc,1,-1 8 if(si(ic).lt.N)goto 9 goto 3 9 do 102 i=ic+1,Nexc 102 si(i)=si(ic)+1 si(ic)=si(ic)+1 goto 5 3 continue c loop end 22222222222222222222222222222222222222222222222222 fcf3=fcsumr return end c ------------------------------------------------------------ function fcf4(nv,mv,nk,w,u,LEXCL,xkr,xki,N7,N,box, 1m1,fcsumi,OH) c fcf1, with an interpolation for k-coefficients implicit none integer*4 maxb,nk0,M3,maxe parameter (maxb=100,M3=100,nk0=200,maxe=200) integer*4 nv(*),mv(*),N7,i,nk,N,Nexc,iex,m1,j, 1state(M3),si(maxe),kstate(M3),ic,LEXCL,km real*8 w(M3,M3),u(*),fcf4,fcsumr,fcsumi,xkr(maxb,2*nk0+1), 1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp, 1ptr,pti,OH,t,xkrn,xkin,ak,dk,dko c fcf Frank-Condon factor c mv,nv excited and ground quantum numbers (1=ground state,nosc=0) c nk number of k-expansion points (k=-nk ... nk) c qs=wql+u c LEXCL maximum sum of k-numbers c xkr,xki real and imaginary parts of the exp. coefficients c N7..N number of oscillators, typically N7=1,N=3*nat-6 c box=L the NM interval (-L/2...L/2) c m1 maximum excitation on one center pi=4.0d0*atan(1.0d0) tp=2.0d0*pi fcsumr=0.0d0 fcsumi=0.0d0 c c distribute LEXCL onto the excited k-expansions N7...N c loop start 22222222222222222222222222222222222222222222222222 do 3 Nexc=0,LEXCL c initial indices - all to the first center: do 41 iex=1,Nexc 41 si(iex)=N7 5 DO 1966 I=N7,N 1966 STATE(I)=0 do 4 iex=1,Nexc 4 state(si(iex))=state(si(iex))+1 do 42 iex=1,Nexc 42 if(state(si(iex)).gt.m1)goto 333 c c odd excitation - positive (N+1)/2 c even excitation - neagtive N/2 do 43 I=N7,N if(mod(state(I),2).eq.0)then kstate(I)=-state(I)/2 else kstate(I)=(state(I)+1)/2 endif if(abs(kstate(i)).gt.nk)goto 333 43 continue c c calculate the phase and product exp(if)xs... fim=0.0d0 do 10 I=N7,N 10 fim=fim+u(I)*kstate(I) fim=fim*tp/box prs=cos(fim) pis=sin(fim) do 11 i=N7,N t = prs*xkr(mv(i),nk+1+kstate(i))-pis*xki(mv(i),nk+1+kstate(i)) pis= prs*xki(mv(i),nk+1+kstate(i))+pis*xkr(mv(i),nk+1+kstate(i)) prs=t 11 if(abs(prs)+abs(pis).lt.OH)goto 333 c c calculate the phase and product (exp(-if)xs*(1)xs*(2)...)xl(1)... ptr=1.0d0 pti=0.0d0 do 111 i=N7,N ak=0.0d0 do 12 j=N7,N c w-first index (j) excited 12 ak=ak+dble(kstate(j))*w(j,i) if(ak.lt.0)then c nearest lower integer: km=int(ak)-1 else km=int(ak) endif if(abs(km).gt.nk)goto 333 if(abs(km+1).gt.nk)goto 333 dk=ak-dble(km) dko=1.0d0-dk c c Fourier coefficients by interpolation: xkrn=xkr(nv(i),nk+1+km)*dko+dk*xkr(nv(i),nk+1+km+1) xkin=xki(nv(i),nk+1+km)*dko+dk*xki(nv(i),nk+1+km+1) t = ptr*xkrn-pti*xkin pti= ptr*xkin+pti*xkrn ptr=t 111 if(abs(ptr)+abs(pti).lt.OH)goto 333 c a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb): fcsumr=fcsumr+(prs*ptr+pis*pti) fcsumi=fcsumi+(prs*pti-pis*ptr) c find index to be changed 333 if(Nexc.eq.0)goto 3 do 8 ic=Nexc,1,-1 8 if(si(ic).lt.N)goto 9 goto 3 9 do 102 i=ic+1,Nexc 102 si(i)=si(ic)+1 si(ic)=si(ic)+1 goto 5 3 continue c loop end 22222222222222222222222222222222222222222222222222 fcf4=fcsumr return end c ------------------------------------------------------------ function fcf5(nv,mv,nk,w,u,LEXCL,xkr,xki,N7,N,box, 1m1,fcsumi,OH) c fcf3 with interpolation for k-coefficients implicit none integer*4 maxb,nk0,M3,maxe parameter (maxb=100,M3=100,nk0=200,maxe=200) integer*4 nv(*),mv(*),N7,i,nk,N,Nexc,iex,m1,j,km, 1state(M3),si(maxe),kstate(M3),ic,LEXCL real*8 w(M3,M3),u(*),fcf5,fcsumr,fcsumi,xkr(maxb,2*nk0+1), 1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp, 1ptr,pti,OH,t,xkrn,xkin,ak,dk,dko c fcf Frank-Condon factor c mv,nv excited and ground quantum numbers (1=ground state,nosc=0) c nk number of k-expansion points (k=-nk ... nk) c qs=wql+u c LEXCL maximum sum of k-numbers c xkr,xki real and imaginary parts of the exp. coefficients c N7..N number of oscillators, typically N7=1,N=3*nat-6 c box=L the NM interval (-L/2...L/2) c m1 maximum excitation on one center pi=4.0d0*atan(1.0d0) tp=2.0d0*pi fcsumr=0.0d0 fcsumi=0.0d0 c distribute k-"excitations" on electronic excited coeff c start with ground as it may need less k-points do 3 Nexc=0,LEXCL c initial indices - all to the first center: do 41 iex=1,Nexc 41 si(iex)=N7 5 DO 1966 I=N7,N 1966 STATE(I)=0 do 4 iex=1,Nexc 4 state(si(iex))=state(si(iex))+1 do 42 iex=1,Nexc 42 if(state(si(iex)).gt.m1)goto 333 c c odd excitation - positive (N+1)/2 c even excitation - neagtive N/2 do 43 I=N7,N if(mod(state(I),2).eq.0)then kstate(I)=-state(I)/2 else kstate(I)=(state(I)+1)/2 endif if(abs(kstate(i)).gt.nk)goto 333 43 continue c c calculate the product xl... prs=1.0d0 pis=0.0d0 do 11 i=N7,N t = prs*xkr(nv(i),nk+1+kstate(i))-pis*xki(nv(i),nk+1+kstate(i)) pis= prs*xki(nv(i),nk+1+kstate(i))+pis*xkr(nv(i),nk+1+kstate(i)) prs=t 11 if(abs(prs)+abs(pis).lt.OH)goto 333 c c calculate the product xs... ptr=1.0d0 pti=0.0d0 fim=0.0d0 do 111 i=N7,N c calculate corresponding ak: ak=0.0d0 do 12 j=N7,N c w-first index excited 12 ak=ak+dble(kstate(j))*w(i,j) if(ak.lt.0)then c nearest lower integer: km=int(ak)-1 else km=int(ak) endif if(abs(km).gt.nk)goto 333 if(abs(km+1).gt.nk)goto 333 dk=ak-dble(km) dko=1.0d0-dk c calculate the phase: fim=fim+u(i)*ak*tp/box c write(6,*)fim,ak,u(i) c c Fourier coefficients by interpolation: xkrn=xkr(mv(i),nk+1+km)*dko+dk*xkr(mv(i),nk+1+km+1) xkin=xki(mv(i),nk+1+km)*dko+dk*xki(mv(i),nk+1+km+1) c t = ptr*xkrn-pti*xkin pti= ptr*xkin+pti*xkrn ptr=t 111 if(abs(ptr)+abs(pti).lt.OH)goto 333 c read(5,*)ic t = ptr*cos(fim)-pti*sin(fim) pti= ptr*sin(fim)+pti*cos(fim) ptr=t c a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb): fcsumr=fcsumr+(prs*ptr+pis*pti) fcsumi=fcsumi+(prs*pti-pis*ptr) c find index to be changed 333 if(Nexc.eq.0)goto 3 do 8 ic=Nexc,1,-1 8 if(si(ic).lt.N)goto 9 goto 3 9 do 102 i=ic+1,Nexc 102 si(i)=si(ic)+1 si(ic)=si(ic)+1 goto 5 3 continue c loop end 22222222222222222222222222222222222222222222222222 fcf5=fcsumr return end c ------------------------------------------------------------ function fcf6(nv,mv,nk,w,u,xkr,xki,N7,N,box, 1fcsumi,OH,time0) c integration by random number generation implicit none integer*4 maxb,nk0,M3 parameter (maxb=100,M3=100,nk0=200) integer*4 nv(*),mv(*),N7,i,nk,N,j,km, 1kstate(M3),k, 1time0,time,mclock real*8 w(M3,M3),u(*),fcf6,fcsumr,fcsumi,xkr(maxb,2*nk0+1), 1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp, 1ptr,pti,OH,t,xkrn,xkin,ak,dk,dko, 1it,ip,im,amax,amaxi,random,y c fcf Frank-Condon factor c mv,nv excited and ground quantum numbers (1=ground state,nosc=0) c nk number of k-expansion points (k=-nk ... nk) c qs=wql+u c LEXCL maximum sum of k-numbers c xkr,xki real and imaginary parts of the exp. coefficients c N7..N number of oscillators, typically N7=1,N=3*nat-6 c box=L the NM interval (-L/2...L/2) c m1 maximum excitation on one center pi=4.0d0*atan(1.0d0) tp=2.0d0*pi fcsumr=0.0d0 fcsumi=0.0d0 time=mclock() c initialize random function do 2993 i=1,44 ak=random(0) 2993 if(ak.gt.100.0d0)write(6,*)'blbost' amax=1.0d0 do 200 i=N7,N amaxi=0.0d0 do 201 k=-nk,nk if(dabs(xkr(nv(i),nk+1+k)).gt.amaxi)amaxi=dabs(xkr(nv(i),nk+1+k)) 201 if(dabs(xkr(mv(i),nk+1+k)).gt.amaxi)amaxi=dabs(xkr(mv(i),nk+1+k)) 200 amax=amax*amaxi**2 it=0 ip=0 im=0 4000 continue do 43 I=N7,N 43 kstate(I)=int(-nk+int(dble(nk+nk)*random(0))) y=-amax+(amax+amax)*random(0) c c calculate the product xl... prs=1.0d0 pis=0.0d0 do 11 i=N7,N t = prs*xkr(nv(i),nk+1+kstate(i))-pis*xki(nv(i),nk+1+kstate(i)) pis= prs*xki(nv(i),nk+1+kstate(i))+pis*xkr(nv(i),nk+1+kstate(i)) prs=t 11 if(abs(prs)+abs(pis).lt.OH)goto 333 c c calculate the product xs... ptr=1.0d0 pti=0.0d0 fim=0.0d0 do 111 i=N7,N c calculate corresponding ak: ak=0.0d0 do 12 j=N7,N c w-first index excited 12 ak=ak+dble(kstate(j))*w(i,j) if(ak.lt.0)then c nearest lower integer: km=int(ak)-1 else km=int(ak) endif if(abs(km).gt.nk)goto 333 if(abs(km+1).gt.nk)goto 333 dk=ak-dble(km) dko=1.0d0-dk c calculate the phase: fim=fim+u(i)*ak*tp/box c c Fourier coefficients by interpolation: xkrn=xkr(mv(i),nk+1+km)*dko+dk*xkr(mv(i),nk+1+km+1) xkin=xki(mv(i),nk+1+km)*dko+dk*xki(mv(i),nk+1+km+1) c t = ptr*xkrn-pti*xkin pti= ptr*xkin+pti*xkrn ptr=t 111 if(abs(ptr)+abs(pti).lt.OH)goto 333 t = ptr*cos(fim)-pti*sin(fim) pti= ptr*sin(fim)+pti*cos(fim) ptr=t c a*.b = (ReaReb+ImaImb) + i (ReaImb-ImaReb): fcsumr=prs*ptr+pis*pti if(y.lt.fcsumr.and.y.gt.0)ip=ip+1.0d0 if(y.gt.fcsumr.and.y.lt.0)im=im+1.0d0 it=it+1.0d0 333 if(mclock()-time.gt.time0)then fcf6=(ip-im)/it*(amax+amax)*dble(nk)**(2*(N-N7+1)) return else goto 4000 endif return end c ------------------------------------------------------------ subroutine fcf7(nv,mv,nk,w,u,LEXCL,xkr,xki,N7,N,box, 1m1,fcsumr,OH,mm) c fcf3 with interpolation for k-coefficients c and with batch of m mm vectors implicit none integer*4 maxb,nk0,M3,maxe,NB0 parameter (maxb=100,M3=100,nk0=200,maxe=200,NB0=200) integer*4 nv(*),mv(*),N7,i,nk,N,Nexc,iex,m1,j,km, 1state(M3),si(maxe),kstate(M3),ic,LEXCL,mm,ii,im real*8 w(M3,M3),u(*),fcsumr(*),xkr(maxb,2*nk0+1), 1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp, 1ptr(NB0),pti(NB0),OH,t,xkrn,xkin,ak,dk,dko c fcf Frank-Condon factor c mv,nv excited and ground quantum numbers (1=ground state,nosc=0) c nk number of k-expansion points (k=-nk ... nk) c qs=wql+u c LEXCL maximum sum of k-numbers c xkr,xki real and imaginary parts of the exp. coefficients c N7..N number of oscillators, typically N7=1,N=3*nat-6 c box=L the NM interval (-L/2...L/2) c m1 maximum excitation on one center pi=4.0d0*atan(1.0d0) tp=2.0d0*pi do 1 i=1,mm 1 fcsumr(i)=0.0d0 c distribute k-"excitations" on electronic excited coeff c start with ground as it may need less k-points do 3 Nexc=0,LEXCL c initial indices - all to the first center: do 41 iex=1,Nexc 41 si(iex)=N7 5 DO 1966 I=N7,N 1966 STATE(I)=0 do 4 iex=1,Nexc 4 state(si(iex))=state(si(iex))+1 do 42 iex=1,Nexc 42 if(state(si(iex)).gt.m1)goto 333 c c odd excitation - positive (N+1)/2 c even excitation - neagtive N/2 do 43 I=N7,N if(mod(state(I),2).eq.0)then kstate(I)=-state(I)/2 else kstate(I)=(state(I)+1)/2 endif if(abs(kstate(i)).gt.nk)goto 333 43 continue c c calculate the product xl...ground state the same for all batch prs=1.0d0 pis=0.0d0 do 11 i=N7,N t = prs*xkr(nv(i),nk+1+kstate(i))-pis*xki(nv(i),nk+1+kstate(i)) pis= prs*xki(nv(i),nk+1+kstate(i))+pis*xkr(nv(i),nk+1+kstate(i)) prs=t 11 if(abs(prs)+abs(pis).lt.OH)goto 333 c c calculate the product xs... do 7 ii=1,mm ptr(ii)=1.0d0 7 pti(ii)=0.0d0 fim=0.0d0 do 111 i=N7,N c calculate corresponding ak: ak=0.0d0 do 12 j=N7,N c w-first index excited 12 ak=ak+dble(kstate(j))*w(i,j) c contribution to the phase: fim=fim+u(i)*ak*tp/box if(ak.lt.0)then c nearest lower integer: km=int(ak)-1 else km=int(ak) endif if(abs(km).gt.nk)goto 333 if(abs(km+1).gt.nk)goto 333 dk=ak-dble(km) dko=1.0d0-dk do 111 ii=1,mm im=mv(ii+(i-1)*NB0) c coefficients by interpolation: xkrn=xkr(im,nk+1+km)*dko+dk*xkr(im,nk+1+km+1) xkin=xki(im,nk+1+km)*dko+dk*xki(im,nk+1+km+1) t = ptr(ii)*xkrn-pti(ii)*xkin pti(ii)= ptr(ii)*xkin+pti(ii)*xkrn 111 ptr(ii)= t do 6 ii=1,mm t = ptr(ii)*cos(fim)-pti(ii)*sin(fim) pti(ii)= ptr(ii)*sin(fim)+pti(ii)*cos(fim) ptr(ii)= t 6 fcsumr(ii)=fcsumr(ii)+(prs*ptr(ii)+pis*pti(ii)) c find index to be changed 333 if(Nexc.eq.0)goto 3 do 8 ic=Nexc,1,-1 8 if(si(ic).lt.N)goto 9 goto 3 9 do 102 i=ic+1,Nexc 102 si(i)=si(ic)+1 si(ic)=si(ic)+1 goto 5 3 continue c loop end 22222222222222222222222222222222222222222222222222 return end c ------------------------------------------------------------ subroutine fcf8(nv,mv,nk,w,u,xkr,xki,N7,N,box, 1fcsumr,OH,time0,mm) c integration by random number generation c and batch of m mm vectors implicit none integer*4 maxb,nk0,M3,NB0 parameter (maxb=100,M3=100,nk0=200,NB0=200) integer*4 nv(*),mv(*),N7,i,nk,N,j,km, 1kstate(M3),ii,k,imm, 1time0,time,mclock,mm real*8 w(M3,M3),u(*),fcsumr(*),sumr, 1xkr(maxb,2*nk0+1), 1xki(maxb,2*nk0+1),prs,pis,box,fim,pi,tp, 1ptr(NB0),pti(NB0),OH,t,xkrn,xkin,ak,dk,dko, 1it(NB0),ip(NB0),im(NB0),amax(NB0),amaxi,random,y(NB0) c fcf Frank-Condon factor c mv,nv excited and ground quantum numbers (1=ground state,nosc=0) c nk number of k-expansion points (k=-nk ... nk) c qs=wql+u c LEXCL maximum sum of k-numbers c xkr,xki real and imaginary parts of the exp. coefficients c N7..N number of oscillators, typically N7=1,N=3*nat-6 c box=L the NM interval (-L/2...L/2) c m1 maximum excitation on one center pi=4.0d0*atan(1.0d0) tp=2.0d0*pi do 1 i=1,mm 1 fcsumr(i)=0.0d0 time=mclock() do 200 ii=1,mm amax(ii)=1.0d0 do 200 i=N7,N amaxi=0.0d0 imm=mv(ii+(i-1)*NB0) do 201 k=-nk,nk if(dabs(xkr(nv(i),nk+1+k)).gt.amaxi)amaxi=dabs(xkr(nv(i),nk+1+k)) 201 if(dabs(xkr(imm,nk+1+k)).gt.amaxi)amaxi=dabs(xkr(imm,nk+1+k)) 200 amax(ii)=amax(ii)*amaxi**2 do 6 ii=1,mm it(ii)=0.0d0 ip(ii)=0.0d0 6 im(ii)=0.0d0 4000 continue do 43 I=N7,N 43 kstate(I)=int(-nk+int(dble(nk+nk)*random(0))) do 8 ii=1,mm 8 y(ii)=-amax(ii)+(amax(ii)+amax(ii))*random(0) c c calculate the product xl... prs=1.0d0 pis=0.0d0 do 11 i=N7,N t = prs*xkr(nv(i),nk+1+kstate(i))-pis*xki(nv(i),nk+1+kstate(i)) pis= prs*xki(nv(i),nk+1+kstate(i))+pis*xkr(nv(i),nk+1+kstate(i)) prs=t 11 if(abs(prs)+abs(pis).lt.OH)goto 333 c c calculate the product xs... do 2 ii=1,mm ptr(ii)=1.0d0 2 pti(ii)=0.0d0 fim=0.0d0 do 111 i=N7,N c calculate corresponding ak: ak=0.0d0 do 12 j=N7,N c w-first index excited 12 ak=ak+dble(kstate(j))*w(i,j) if(ak.lt.0)then c nearest lower integer: km=int(ak)-1 else km=int(ak) endif if(abs(km).gt.nk)goto 333 if(abs(km+1).gt.nk)goto 333 dk=ak-dble(km) dko=1.0d0-dk c calculate the phase: fim=fim+u(i)*ak*tp/box c do 111 ii=1,mm imm=mv(ii+(i-1)*NB0) xkrn=xkr(imm,nk+1+km)*dko+dk*xkr(imm,nk+1+km+1) xkin=xki(imm,nk+1+km)*dko+dk*xki(imm,nk+1+km+1) t =ptr(ii)*xkrn-pti(ii)*xkin pti(ii)=ptr(ii)*xkin+pti(ii)*xkrn 111 ptr(ii)=t do 4 ii=1,mm t = ptr(ii)*cos(fim)-pti(ii)*sin(fim) pti(ii)= ptr(ii)*sin(fim)+pti(ii)*cos(fim) ptr(ii)= t sumr=prs*ptr(ii)+pis*pti(ii) if(y(ii).lt.sumr.and.y(ii).gt.0)ip(ii)=ip(ii)+1.0d0 if(y(ii).gt.sumr.and.y(ii).lt.0)im(ii)=im(ii)+1.0d0 4 it(ii)=it(ii)+1.0d0 333 if(mclock()-time.gt.time0)then do 5 ii=1,mm 5 fcsumr(ii)= 1 (ip(ii)-im(ii))/it(ii)*amax(ii)*2.0d0*dble(nk)**(2*(N-N7+1)) return else goto 4000 endif return end c ------------------------------------------------------------ c ------------------------------------------------------------ subroutine intq(n,k,r,i) c int(q^n exp(-q^2/2-ik.q)/sqrt(2*pi),q=-inf...inf) c r,i .. real and imag. parts implicit none integer*4 n,ni real*8 k,r,i,anir,anii,anipr,anipi,anippr,anippi ni=0 anir=1.0d0 anii=0.0d0 if(n.eq.ni)then r=anir*exp(-k**2/2.0d0) i=anii return endif anipr=0.0d0 anipi=-k if(n.eq.ni+1)then r=anipr i=anipi*exp(-k**2/2.0d0) return endif c I(n+2) = -i k I(n+1) + (n+1) I(n): 1 anippr= k*anipi+dble(ni+1)*anir anippi=-k*anipr+dble(ni+1)*anii if(n.eq.ni+2)then r=anippr*exp(-k**2/2.0d0) i=anippi*exp(-k**2/2.0d0) return else ni=ni+1 anir=anipr anii=anipi anipr=anippr anipi=anippi goto 1 endif end c ------------------------------------------------------------ subroutine report(s) character*(*) s write(6,*)s stop end c ------------------------------------------------------------ function fik(k,x) implicit none real*8 fik,x,w,hk,p,ank,ff integer maxb,M3 parameter (maxb=100,M3=100) common/ho/p(maxb,maxb),ank(maxb),w(M3) integer*4 k,kk ff=ank(k)*exp(-x**2/2.0d0) hk=0.0d0 do 1 kk=1,k 1 hk=hk+p(k,kk)*x**(kk-1) fik=ff*hk return end c ------------------------------------------------------------ c ====================================================================== subroutine inv(fx,M,M0,ipiv,work,a) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION fx(M0,M0),a(M0,M0),ipiv(M0),work(3*M0*M0) do 3 i=1,M do 3 j=i,M a(j,i)=fx(i,j) 3 a(i,j)=fx(i,j) M1=M M2=M call dgetrf(M1,M2,a,M0,ipiv,info) if(info.eq.0)then write(6,*)' Pivoting OK' else call report('pivoting failed') endif c lwork=3*M*M call dgetri(M,a,M0,ipiv,work,lwork,info) c if(info.eq.0)write(6,*)'successful inversion' if(info.lt.0)then call report('invalid argument') endif c tol=1.0d-7 nmo=M do i=1,nmo do j=1,nmo sum=0.0d0 sum2=0.0d0 do ii=1,nmo sum=sum+a(i,ii)*fx(ii,j) sum2=sum2+a(ii,i)*fx(j,ii) enddo ic=0 if(abs(sum ).gt.tol.and.i.ne.j)ic=1 if(abs(sum2).gt.tol.and.i.ne.j)ic=2 if(abs(sum -1.0d0).gt.tol.and.i.eq.j)ic=3 if(abs(sum2-1.0d0).gt.tol.and.i.eq.j)ic=4 if(ic.ne.0)then write(6,6990)((a(i1,j1),j1=1,M),i1=1,M) write(6,*)'ic = ',ic write(6,*)i,j,sum,sum2 write(6,6990)((fx(i1,j1),j1=1,M),i1=1,M) 6990 format(7f10.4) call report('Inversion error') endif enddo enddo c write(6,*)'Inversion Check OK' return end c ================================== subroutine ALLIGN(N7,N,W,UU,e1,e2,Wi) implicit none INTEGER*4 NTA,MENDELEV,MXB,i,ia,ix,iy,iz,ii,NAT, 1ierr,iter,j,MB,maxiter,itrial, 2IX3,M3,NM,iq,idum,iix,N7,N PARAMETER (M3=100,NTA=M3/3,MENDELEV=89,MXB=25,maxiter=1000) real*8 pi,A,XS,XB,RTOL,am, 1oxx,amxx,ck,sk,dr1,dr2,fi,fk,pp,rat,wx,yy,z0, 2tol,rr(NTA,3),c,s2r(M3,M3),jmat(M3,M3),kmat(M3),Wi(M3,M3), 3X(NTA),r(NTA,3),rg(NTA,3),d(NTA,3),u(3),re(NTA,3) integer*4 ITY(NTA) real*4 amas dimension amas(MENDELEV) data amas/1.007825,4.003, 2 6.941, 9.012, 10.810,12.0 ,14.003074,15.994915,18.9984,20.179, 3 22.990,24.305, 26.981,28.086,30.9737,31.972,34.968852,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ data z0/0.0d0/ COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,MB,ITER real*8 s1(M3,M3),s2(M3,M3),e1(M3),e2(M3),UU(*),W(M3,M3), 1tp,clight,const,hpl,nav logical lex C pi=4.0d0*atan(1.0d0) tol=1.0d-7 tp=2.0d0*pi OPEN(4,FILE='fg.inp',STATUS='OLD') READ(4,*)NM,idum,NAT if(NAT.gt.NTA)call report('too many atoms') if(NM.ne.3*NAT)call report('3xNat x 3xNat S matrix required') DO 7 I=1,NAT READ(4,*)ITY(I),(rg(I,ix),ix=1,3) 7 X(I)=dble(AMAS(ITY(I))) read(4,*) do 1 I=1,NAT do 1 iq=1,NM 1 read(4,*)idum,idum,(s1(3*(I-1)+ix,iq),ix=1,3) read(4,*) read(4,*)(e1(iq),iq=1,NM) close(4) CLOSE(4) inquire(file='FTRY.INP',exist=lex) if(lex)then open(4,file='FTRY.INP') read(4,*) read(4,*) read(4,444)(X(I),I=1,NAT) read(4,*) read(4,444)(X(I),I=1,NAT) 444 format(6F12.6) close(4) write(6,*)'Masses read from FTRY.INP' endif c AM=z0 do 3 i=1,nat 3 AM=AM+X(i) WRITE(6,609)NAT,AM 609 format(' fg.inp read,',i4,' ATOMS, MOLECULAR MASS ',f10.3) OPEN(4,FILE='fe.inp',STATUS='OLD') READ(4,*) DO 72 I=1,NAT 72 READ(4,*)ITY(I),(re(I,ix),ix=1,3) read(4,*) do 2 I=1,NAT do 2 iq=1,NM 2 read(4,*)idum,idum,(s2(3*(I-1)+ix,iq),ix=1,3) read(4,*) read(4,*)(e2(iq),iq=1,NM) close(4) CLOSE(4) write(6,*)'fe.inp read' do 303 ix=1,3 c=z0 do 300 ia=1,NAT 300 c = c+rg(ia,ix)*X(ia) do 303 ia=1,NAT 303 rg(ia,ix)=rg(ia,ix)-c/AM write(6,*)'fg put to origin' do 203 ix=1,3 c=z0 do 200 ia=1,NAT 200 c=c+re(ia,ix)*X(ia) do 203 ia=1,NAT 203 re(ia,ix)=re(ia,ix)-c/AM write(6,*)'fe put to origin' c c rotate re to rg do 2031 ia=1,NAT do 2031 ix=1,3 2031 r(ia,ix)=re(ia,ix) itrial=1 88881 iter=0 8888 iter=iter+1 if(iter.gt.maxiter)then write(6,*)'Maximum number of iterations exceeded' itrial=itrial+1 if(itrial.lt.4)then write(6,*)'Trial ',itrial goto 88881 endif write(6,*)'Iteration did not converge' goto 881 endif do 210 ia=1,NAT do 210 ix=1,3 210 d(ia,ix)=r(ia,ix)-rg(ia,ix) c c loop over the three axes,(try three different orderings): do 208 IX3=1,3 IX=itrial+IX3-1 if(IX.eq.4)IX=1 if(IX.eq.5)IX=2 IY=IX+1 if(IY.eq.4)IY=1 IZ=IY+1 if(IZ.eq.4)IZ=1 oxx=z0 amxx=z0 do 209 ia=1,NAT dr1=dsqrt(rg(ia,IY)**2+rg(ia,IZ)**2) dr2=dsqrt( r(ia,IY)**2+ r(ia,IZ)**2) pp=rg(ia,IY)*r(ia,IY)+rg(ia,IZ)*r(ia,IZ) if(dabs(dr1*dr2).gt.1.0d-8)then rat=pp/dr1/dr2 if(rat.gt. 1.0d0-1.0d-13)then fi=z0 else if(rat.lt.-1.0d0+1.0d-13)then fi=pi else fi=acos(rat) endif endif wx=rg(ia,IY)*d(ia,IZ)-rg(ia,IZ)*d(ia,IY) if(wx.lt.0.0d0)fi=-fi c amxx=amxx+X(ia)*fi*(5.0d0*dr1**2+2.0d0*dr2**2-dr1*dr2)/6.0d0 amxx=amxx+X(ia)*fi*(dr1**2+dr2**2+dr1*dr2)/3.0d0 oxx=oxx+X(ia)*dr2**2 endif 209 continue if(dabs(oxx).gt.1.0d-8)then fk=amxx/oxx else fk=z0 endif u(IX)=fk c c rotate back by this angle: ck=cos(-fk) sk=sin(-fk) do 208 ia=1,NAT YY=r(ia,IY) r(ia,IY)=YY*ck-r(ia,IZ)*sk 208 r(ia,IZ)=YY*sk+r(ia,IZ)*ck c write(6,*)iter,u(1)**2+u(2)**2+u(3)**2 if(u(1)**2+u(2)**2+u(3)**2.gt.tol)goto 8888 881 write(6,608)iter 608 format(i7,' iterations') open(41,file='tr.x') write(41,*)'rotated coordinates of exc' write(41,*)NAT do 2043 ia=1,NAT 2043 write(41,4141)ITY(ia),(r(ia,ix),ix=1,3) close(41) write(6,*) MB=NAT do 217 ii=1,MB do 217 ix=1,3 XB(ii,ix)=r(ii,ix) 217 XS(ii,ix)=re(ii,ix) call DOMATRIX(IERR) c XB = A XS c rotated Xe = A original Xe write(6,6004)A 6004 format(' Transformation matrix:',/,3(3F10.4,/)) c mass-weighed S-matrices do 214 ia=1,NAT do 214 ix=1,3 ii=3*(ia-1)+ix do 214 iq=1,NM s1(ii,iq)=s1(ii,iq)*dsqrt(X(ia)) 214 s2(ii,iq)=s2(ii,iq)*dsqrt(X(ia)) c now transpose is inverted do 213 ia=1,NAT do 213 ix=1,3 ii=3*(ia-1)+ix rr(ia,ix)=0.0d0 do 2131 iix=1,3 2131 rr(ia,ix)=rr(ia,ix)+A(ix,iix)*r(ia,iix) do 213 iq=1,NM s2r(ii,iq)=0.0d0 do 213 iix=1,3 213 s2r(ii,iq)=s2r(ii,iq)+A(ix,iix)*s2(3*(ia-1)+iix,iq) write(6,*)'fe transformed' open(41,file='xr.x') write(41,*)'rotated coordinates of exc' write(41,*)NAT do 2042 ia=1,NAT 2042 write(41,4141)ITY(ia),(rr(ia,ix),ix=1,3) 4141 format(i4,3f12.6) close(41) open(41,file='xg.x') write(41,*)'ground geometry in origin' write(41,*)NAT do 204 ia=1,NAT 204 write(41,4141)ITY(ia),(rg(ia,ix),ix=1,3) close(41) do 215 i=1,NM do 215 j=1,NM w(i,j)=0.0d0 wi(i,j)=0.0d0 jmat(i,j)=0.0d0 do 215 ii=1,3*nat c w,wi - first index always excited c transformation matrix between dimensionless normal mode coords jmat(i,j)=jmat(i,j)+s2r(ii,i)*s1(ii,j) wi(i,j)=wi(i,j)+s2r(ii,i)*dsqrt(e1(j)/e2(i))*s1(ii,j) 215 w(i,j)=w(i,j)+s2r(ii,i)*dsqrt(e2(i)/e1(j))*s1(ii,j) clight=3.0d10 nav=6.022d23 hpl=1.054d-34 const=dsqrt(1.0d-3*tp*clight/nav/hpl)*1.0d-10 do 218 i=1,NM uu(i)=0.0d0 kmat(i)=0.0d0 do 218 ia=1,nat do 218 ix=1,3 ii=3*(ia-1)+ix kmat(i)=kmat(i)-s2r(ii,i)*(rg(ia,ix)-rr(ia,ix))*dsqrt(X(ia)) 218 uu(i)=uu(i)+ 1dsqrt(e2(i)*X(ia))*s2r(ii,i)*const*(rg(ia,ix)-rr(ia,ix)) write(6,6007) 6007 format(/,' J-matrix:') do 221 j=1,3*NAT-6 221 write(6,6008)(jmat(i,j),i=1,3*NAT-6) write(6,6009) 6009 format(/,' K-vector:') write(6,6008)(kmat(i),i=1,3*NAT-6) 6008 format(10f8.2) write(6,6005) 6005 format(' mode Eg Ee shift',/,'----------------') NM=3*NAT-6 do 219 i=1,NM 219 write(6,6006)i,e1(i),e2(i),uu(i) 6006 format(i5,2f8.1,f9.3) N7=1 N=NM return END SUBROUTINE DOMATRIX(IERR) C This is a Fortran version of Petr Malon's subroutine Amoeba IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=25) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER IERR=0 C STARTING VALUES FOR THE ITERATION: ANG(1)=0.2d0 ANG(2)=0.0d0 ANG(3)=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(ANG) ANG(2)=ANG(1) ANG(1)=0.0d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(ANG) ANG(3)=ANG(2) ANG(2)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(ANG) ANG(1)=ANG(3) ANG(2)=ANG(3) ANG(3)=0.0d0 DO 4 I=1,3 4 P(4,I)=ANG(I) Y(4)=FU(ANG) FTOL=1.0d-7 ITMAX=1500 ITER=0 99999 ILO=1 IF(Y(1).GT.Y(2))THEN IHI=1 INHI=2 ELSE IHI=2 INHI=1 ENDIF DO 5 I=1,4 IF(Y(I).LT.Y(ILO))ILO=I IF(Y(I).GT.Y(IHI))THEN INHI=IHI IHI=I ELSE IF(Y(I).GT.Y(INHI))THEN IF(I.NE.IHI)INHI=I ENDIF ENDIF 5 CONTINUE IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.1.0d-11)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)RETURN IF(ITER.EQ.ITMAX)THEN WRITE(3,*)' Rotation has not converged !' IERR=1 RETURN ENDIF ITER=ITER+1 DO 55 I=1,3 55 PBAR(I)=0.0d0 DO 6 I=1,4 IF(I.NE.IHI)THEN DO 7 J=1,3 7 PBAR(J)=PBAR(J)+P(I,J) ENDIF 6 CONTINUE DO 8 J=1,3 PBAR(J)=PBAR(J)/3.0d0 8 PR(J)=2.0d0*PBAR(J)-P(IHI,J) YPR=FU(PR) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,3 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(ILO))THEN DO 10 J=1,3 10 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 11 J=1,3 11 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ELSE IF(YPR.GE.Y(INHI))THEN IF(YPR.LT.Y(IHI))THEN DO 12 J=1,3 12 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF DO 13 J=1,3 13 PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(IHI))THEN DO 14 J=1,3 14 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 15 I=1,4 IF(I.NE.ILO)THEN DO 16 J=1,3 PR(J)=0.5d0*(P(I,J)+P(ILO,J)) 16 P(I,J)=PR(J) Y(I)=FU(PR) ENDIF 15 CONTINUE ENDIF ELSE DO 17 J=1,3 17 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ENDIF GOTO 99999 END FUNCTION FU(ANG) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=25) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER DIMENSION ANG(3) S1=SIN(ANG(1)) S2=SIN(ANG(2)) S3=SIN(ANG(3)) C1=COS(ANG(1)) C2=COS(ANG(2)) C3=COS(ANG(3)) C ang1 = theta C ang2 = phi C ang3 = cappa A(1,1)=C1*C2*C3-S2*S3 A(1,2)=C1*S2*C3+C2*S3 A(1,3)=-S1*C3 A(2,1)=-C1*C2*S3-S2*C3 A(2,2)=-C1*S2*S3+C2*C3 A(2,3)=S1*S3 A(3,1)=S1*C2 A(3,2)=S1*S2 A(3,3)=C1 R=0.0d0 DO 1 I=1,NAT TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3) TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3) TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3) 1 R=R+(TX-XB(I,1))**2+(TY-XB(I,2))**2+(TZ-XB(I,3))**2 FU=R RETURN END c c =========================================== c LAPACK ROUTINES: c =========================================== SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1999 * * .. Scalar Arguments .. INTEGER*4 INFO, LDA, LWORK, N * .. * .. Array Arguments .. INTEGER*4 IPIV( * ) DOUBLE PRECISION A( LDA, * ), WORK( * ) * .. * * Purpose * ======= * * DGETRI computes the inverse of a matrix using the LU factorization * computed by DGETRF. * * This method inverts U and then computes inv(A) by solving the system * inv(A)*L = inv(U) for inv(A). * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the factors L and U from the factorization * A = P*L*U as computed by DGETRF. * On exit, if INFO = 0, the inverse of the original matrix A. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * IPIV (input) INTEGER*4 array, dimension (N) * The pivot indices from DGETRF; for 1<=i<=N, row i of the * matrix was interchanged with row IPIV(i). * * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) * On exit, if INFO=0, then WORK(1) returns the optimal LWORK. * * LWORK (input) INTEGER * The dimension of the array WORK. LWORK >= max(1,N). * For optimal performance LWORK >= N*NB, where NB is * the optimal blocksize returned by ILAENV. * * If LWORK = -1, then a workspace query is assumed; the routine * only calculates the optimal size of the WORK array, returns * this value as the first entry of the WORK array, and no error * message related to LWORK is issued by XERBLA. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, U(i,i) is exactly zero; the matrix is * singular and its inverse could not be computed. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL LQUERY INTEGER*4 I, IWS, J, JB, JJ, JP, LDWORK, LWKOPT, NB, $ NBMIN, NN * .. * .. External Functions .. INTEGER*4 ILAENV EXTERNAL ILAENV * .. * .. External Subroutines .. EXTERNAL DGEMM, DGEMV, DSWAP, DTRSM, DTRTRI, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 NB = ILAENV( 1, 'DGETRI', ' ', N, -1, -1, -1 ) LWKOPT = N*NB WORK( 1 ) = LWKOPT LQUERY = ( LWORK.EQ.-1 ) IF( N.LT.0 ) THEN INFO = -1 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -3 ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN INFO = -6 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETRI', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * * Form inv(U). If INFO > 0 from DTRTRI, then U is singular, * and the inverse is not computed. * CALL DTRTRI( 'Upper', 'N', N, A, LDA, INFO ) IF( INFO.GT.0 ) $ RETURN * NBMIN = 2 LDWORK = N IF( NB.GT.1 .AND. NB.LT.N ) THEN IWS = MAX( LDWORK*NB, 1 ) IF( LWORK.LT.IWS ) THEN NB = LWORK / LDWORK NBMIN = MAX( 2, ILAENV( 2, 'DGETRI', ' ', N, -1, -1, -1 ) ) END IF ELSE IWS = N END IF * * Solve the equation inv(A)*L = inv(U) for inv(A). * IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN * * Use unblocked code. * DO 20 J = N, 1, -1 * * Copy current column of L to WORK and replace with zeros. * DO 10 I = J + 1, N WORK( I ) = A( I, J ) A( I, J ) = ZERO 10 CONTINUE * * Compute current column of inv(A). * IF( J.LT.N ) $ CALL DGEMV( 'N', N, N-J, -ONE, A( 1, J+1 ), $ LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 ) 20 CONTINUE ELSE * * Use blocked code. * NN = ( ( N-1 ) / NB )*NB + 1 DO 50 J = NN, 1, -NB JB = MIN( NB, N-J+1 ) * * Copy current block column of L to WORK and replace with * zeros. * DO 40 JJ = J, J + JB - 1 DO 30 I = JJ + 1, N WORK( I+( JJ-J )*LDWORK ) = A( I, JJ ) A( I, JJ ) = ZERO 30 CONTINUE 40 CONTINUE * * Compute current block column of inv(A). * IF( J+JB.LE.N ) $ CALL DGEMM( 'N', 'N', N, JB, $ N-J-JB+1, -ONE, A( 1, J+JB ), LDA, $ WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA ) CALL DTRSM( 'R', 'L', 'N', 'U', N, JB, $ ONE, WORK( J ), LDWORK, A( 1, J ), LDA ) 50 CONTINUE END IF * * Apply column interchanges. * DO 60 J = N - 1, 1, -1 JP = IPIV( J ) IF( JP.NE.J ) $ CALL DSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 ) 60 CONTINUE * WORK( 1 ) = IWS RETURN * * End of DGETRI * END SUBROUTINE DTRTI2( UPLO, DIAG, N, A, LDA, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER DIAG, UPLO INTEGER*4 INFO, LDA, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DTRTI2 computes the inverse of arreal upper or lower triangular * matrix. * * This is the Level 2 BLAS version of the algorithm. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * Specifies whether the matrix A is upper or lower triangular. * = 'U': Upper triangular * = 'L': Lower triangular * * DIAG (input) CHARACTER*1 * Specifies whether or not the matrix A is unit triangular. * = 'N': Non-unit triangular * = 'U': Unit triangular * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the triangular matrix A. If UPLO = 'U', the * leading n by n upper triangular part of the array A contains * the upper triangular matrix, and the strictly lower * triangular part of A is not referenced. If UPLO = 'L', the * leading n by n lower triangular part of the array A contains * the lower triangular matrix, and the strictly upper * triangular part of A is not referenced. If DIAG = 'U', the * diagonal elements of A are also not referenced and are * assumed to be 1. * * On exit, the (triangular) inverse of the original matrix, in * the same storage format. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL NOUNIT, UPPER INTEGER*4 J DOUBLE PRECISION AJJ * .. * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. * .. External Subroutines .. EXTERNAL DSCAL, DTRMV, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LBAME( UPLO, 'U' ) NOUNIT = LBAME( DIAG, 'N' ) IF( .NOT.UPPER .AND. .NOT.LBAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( .NOT.NOUNIT .AND. .NOT.LBAME( DIAG, 'U' ) ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DTRTI2', -INFO ) RETURN END IF * IF( UPPER ) THEN * * Compute inverse of upper triangular matrix. * DO 10 J = 1, N IF( NOUNIT ) THEN A( J, J ) = ONE / A( J, J ) AJJ = -A( J, J ) ELSE AJJ = -ONE END IF * * Compute elements 1:j-1 of j-th column. * CALL DTRMV( 'Upper', 'N', DIAG, J-1, A, LDA, $ A( 1, J ), 1 ) CALL DSCAL( J-1, AJJ, A( 1, J ), 1 ) 10 CONTINUE ELSE * * Compute inverse of lower triangular matrix. * DO 20 J = N, 1, -1 IF( NOUNIT ) THEN A( J, J ) = ONE / A( J, J ) AJJ = -A( J, J ) ELSE AJJ = -ONE END IF IF( J.LT.N ) THEN * * Compute elements j+1:n of j-th column. * CALL DTRMV( 'L', 'N', DIAG, N-J, $ A( J+1, J+1 ), LDA, A( J+1, J ), 1 ) CALL DSCAL( N-J, AJJ, A( J+1, J ), 1 ) END IF 20 CONTINUE END IF * RETURN * * End of DTRTI2 * END SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * March 31, 1993 * * .. Scalar Arguments .. CHARACTER DIAG, UPLO INTEGER*4 INFO, LDA, N * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ) * .. * * Purpose * ======= * * DTRTRI computes the inverse of arreal upper or lower triangular * matrix A. * * This is the Level 3 BLAS version of the algorithm. * * Arguments * ========= * * UPLO (input) CHARACTER*1 * = 'U': A is upper triangular; * = 'L': A is lower triangular. * * DIAG (input) CHARACTER*1 * = 'N': A is non-unit triangular; * = 'U': A is unit triangular. * * N (input) INTEGER * The order of the matrix A. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the triangular matrix A. If UPLO = 'U', the * leading N-by-N upper triangular part of the array A contains * the upper triangular matrix, and the strictly lower * triangular part of A is not referenced. If UPLO = 'L', the * leading N-by-N lower triangular part of the array A contains * the lower triangular matrix, and the strictly upper * triangular part of A is not referenced. If DIAG = 'U', the * diagonal elements of A are also not referenced and are * assumed to be 1. * On exit, the (triangular) inverse of the original matrix, in * the same storage format. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,N). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, A(i,i) is exactly zero. The triangular * matrix is singular and its inverse can not be computed. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. LOGICAL NOUNIT, UPPER INTEGER*4 J, JB, NB, NN * .. * .. External Functions .. LOGICAL LBAME INTEGER*4 ILAENV EXTERNAL LBAME, ILAENV * .. * .. External Subroutines .. EXTERNAL DTRMM, DTRSM, DTRTI2, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LBAME( UPLO, 'U' ) NOUNIT = LBAME( DIAG, 'N' ) IF( .NOT.UPPER .AND. .NOT.LBAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( .NOT.NOUNIT .AND. .NOT.LBAME( DIAG, 'U' ) ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DTRTRI', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 ) $ RETURN * * Check for singularity if non-unit. * IF( NOUNIT ) THEN DO 10 INFO = 1, N IF( A( INFO, INFO ).EQ.ZERO ) $ RETURN 10 CONTINUE INFO = 0 END IF * * Determine the block size for this environment. * NB = ILAENV( 1, 'DTRTRI', UPLO // DIAG, N, -1, -1, -1 ) IF( NB.LE.1 .OR. NB.GE.N ) THEN * * Use unblocked code * CALL DTRTI2( UPLO, DIAG, N, A, LDA, INFO ) ELSE * * Use blocked code * IF( UPPER ) THEN * * Compute inverse of upper triangular matrix * DO 20 J = 1, N, NB JB = MIN( NB, N-J+1 ) * * Compute rows 1:j-1 of current block column * CALL DTRMM( 'L', 'Upper', 'N', DIAG, J-1, $ JB, ONE, A, LDA, A( 1, J ), LDA ) CALL DTRSM( 'R', 'Upper', 'N', DIAG, J-1, $ JB, -ONE, A( J, J ), LDA, A( 1, J ), LDA ) * * Compute inverse of current diagonal block * CALL DTRTI2( 'U', DIAG, JB, A( J, J ), LDA, INFO ) 20 CONTINUE ELSE * * Compute inverse of lower triangular matrix * NN = ( ( N-1 ) / NB )*NB + 1 DO 30 J = NN, 1, -NB JB = MIN( NB, N-J+1 ) IF( J+JB.LE.N ) THEN * * Compute rows j+jb:n of current block column * CALL DTRMM( 'L', 'L', 'N', DIAG, $ N-J-JB+1, JB, ONE, A( J+JB, J+JB ), LDA, $ A( J+JB, J ), LDA ) CALL DTRSM( 'R', 'L', 'N', DIAG, $ N-J-JB+1, JB, -ONE, A( J, J ), LDA, $ A( J+JB, J ), LDA ) END IF * * Compute inverse of current diagonal block * CALL DTRTI2( 'L', DIAG, JB, A( J, J ), LDA, INFO ) 30 CONTINUE END IF END IF * RETURN * * End of DTRTRI * END LOGICAL FUNCTION LBAME( CA, CB ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER CA, CB * .. * * Purpose * ======= * * LBAME returns .TRUE. if CA is the same letter as CB regardless of * case. * * Arguments * ========= * * CA (input) CHARACTER*1 * CB (input) CHARACTER*1 * CA and CB specify the single characters to be compared. * * ===================================================================== * * .. Intrinsic Functions .. INTRINSIC ICHAR * .. * .. Local Scalars .. INTEGER*4 INTA, INTB, ZCODE * .. * .. Executable Statements .. * * Test if the characters are equal * LBAME = CA.EQ.CB IF( LBAME ) $ RETURN * * Now test for equivalence if both characters are alphabetic. * ZCODE = ICHAR( 'Z' ) * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LBAME = INTA.EQ.INTB * * RETURN * * End of LBAME * END SUBROUTINE XERBLA( SRNAME, INFO ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * September 30, 1994 * * .. Scalar Arguments .. CHARACTER*6 SRNAME INTEGER*4 INFO * .. * * Purpose * ======= * * XERBLA is an error handler for the LAPACK routines. * It is called by an LAPACK routine if an input parameter has an * invalid value. A message is printed and execution stops. * * Installers may consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Arguments * ========= * * SRNAME (input) CHARACTER*6 * The name of the routine which called XERBLA. * * INFO (input) INTEGER * The position of the invalid parameter in the parameter list * of the calling routine. * * ===================================================================== * * .. Executable Statements .. * WRITE( *, FMT = 9999 )SRNAME, INFO * STOP * 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ', $ 'an illegal value' ) * * End of XERBLA * END SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 TRANSA, TRANSB INTEGER*4 M, N, K, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * Purpose * ======= * * DGEMM performs one of the matrix-matrix operations * * C := alpha*op( A )*op( B ) + beta*C, * * where op( X ) is one of * * op( X ) = X or op( X ) = X', * * alpha and beta are scalars, and A, B and C are matrices, with op( A ) * an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. * * Parameters * ========== * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n', op( A ) = A. * * TRANSA = 'T' or 't', op( A ) = A'. * * TRANSA = 'C' or 'c', op( A ) = A'. * * Unchanged on exit. * * TRANSB - CHARACTER*1. * On entry, TRANSB specifies the form of op( B ) to be used in * the matrix multiplication as follows: * * TRANSB = 'N' or 'n', op( B ) = B. * * TRANSB = 'T' or 't', op( B ) = B'. * * TRANSB = 'C' or 'c', op( B ) = B'. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix * op( A ) and of the matrix C. M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix * op( B ) and the number of columns of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of columns of the matrix * op( A ) and the number of rows of the matrix op( B ). K must * be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is * k when TRANSA = 'N' or 'n', and is m otherwise. * Before entry with TRANSA = 'N' or 'n', the leading m by k * part of the array A must contain the matrix A, otherwise * the leading k by m part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANSA = 'N' or 'n' then * LDA must be at least max( 1, m ), otherwise LDA must be at * least max( 1, k ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is * n when TRANSB = 'N' or 'n', and is k otherwise. * Before entry with TRANSB = 'N' or 'n', the leading k by n * part of the array B must contain the matrix B, otherwise * the leading n by k part of the array B must contain the * matrix B. * Unchanged on exit. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. When TRANSB = 'N' or 'n' then * LDB must be at least max( 1, k ), otherwise LDB must be at * least max( 1, n ). * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then C need not be set on input. * Unchanged on exit. * * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). * Before entry, the leading m by n part of the array C must * contain the matrix C, except when beta is zero, in which * case C need not be set on entry. * On exit, the array C is overwritten by the m by n matrix * ( alpha*op( A )*op( B ) + beta*C ). * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL NOTA, NOTB INTEGER*4 I, INFO, J, L, NROWA, NROWB DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Set NOTA and NOTB as true if A and B respectively are not * transposed and set NROWA, NCOLA and NROWB as the number of rows * and columns of A and the number of rows of B respectively. * NOTA = LBAME( TRANSA, 'N' ) NOTB = LBAME( TRANSB, 'N' ) IF( NOTA )THEN NROWA = M c NCOLA = K ELSE NROWA = K c NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * Test the input parameters. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.LBAME( TRANSA, 'C' ) ).AND. $ ( .NOT.LBAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.LBAME( TRANSB, 'C' ) ).AND. $ ( .NOT.LBAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And if alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * Start the operations. * IF( NOTB )THEN IF( NOTA )THEN * * Form C := alpha*A*B + beta*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE * * Form C := alpha*A'*B + beta*C * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN * * Form C := alpha*A*B' + beta*C * DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE * * Form C := alpha*A'*B' + beta*C * DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF * RETURN * * End of DGEMM . * END SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA, BETA INTEGER*4 INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * X - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, the incremented array Y * must contain the vector y. On exit, Y is overwritten by the * updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER*4 I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LBAME( TRANS, 'N' ).AND. $ .NOT.LBAME( TRANS, 'T' ).AND. $ .NOT.LBAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Set LENX and LENY, the lengths of the vectors x and y, and set * up the start points in X and Y. * IF( LBAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN IF( LBAME( TRANS, 'N' ) )THEN * * Form y := alpha*A*x + y. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF ELSE * * Form y := alpha*A'*x + y. * JY = KY IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = ZERO DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE DO 120, J = 1, N TEMP = ZERO IX = KX DO 110, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * End of DGEMV . * 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 modified 12/3/93, array(1) declarations changed to array(*) c real*8 dx(*),dy(*),dtemp integer*4 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 SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER*4 M, N, LDA, LDB DOUBLE PRECISION ALPHA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DTRMM performs one of the matrix-matrix operations * * B := alpha*op( A )*B, or B := alpha*B*op( A ), * * where alpha is a scalar, B is an m by n matrix, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) multiplies B from * the left or right as follows: * * SIDE = 'L' or 'l' B := alpha*op( A )*B. * * SIDE = 'R' or 'r' B := alpha*B*op( A ). * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the matrix B, and on exit is overwritten by the * transformed matrix. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER*4 I, INFO, J, K, NROWA DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LBAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LBAME( DIAG , 'N' ) UPPER = LBAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LBAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LBAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LBAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LBAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LBAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LBAME( DIAG , 'U' ) ).AND. $ ( .NOT.LBAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRMM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LBAME( TRANSA, 'N' ) )THEN * * Form B := alpha*A*B. * IF( UPPER )THEN DO 50, J = 1, N DO 40, K = 1, M IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) DO 30, I = 1, K - 1 B( I, J ) = B( I, J ) + TEMP*A( I, K ) 30 CONTINUE IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) B( K, J ) = TEMP END IF 40 CONTINUE 50 CONTINUE ELSE DO 80, J = 1, N DO 70 K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) B( K, J ) = TEMP IF( NOUNIT ) $ B( K, J ) = B( K, J )*A( K, K ) DO 60, I = K + 1, M B( I, J ) = B( I, J ) + TEMP*A( I, K ) 60 CONTINUE END IF 70 CONTINUE 80 CONTINUE END IF ELSE * * Form B := alpha*A'*B. * IF( UPPER )THEN DO 110, J = 1, N DO 100, I = M, 1, -1 TEMP = B( I, J ) IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 90, K = 1, I - 1 TEMP = TEMP + A( K, I )*B( K, J ) 90 CONTINUE B( I, J ) = ALPHA*TEMP 100 CONTINUE 110 CONTINUE ELSE DO 140, J = 1, N DO 130, I = 1, M TEMP = B( I, J ) IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 120, K = I + 1, M TEMP = TEMP + A( K, I )*B( K, J ) 120 CONTINUE B( I, J ) = ALPHA*TEMP 130 CONTINUE 140 CONTINUE END IF END IF ELSE IF( LBAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*A. * IF( UPPER )THEN DO 180, J = N, 1, -1 TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 150, I = 1, M B( I, J ) = TEMP*B( I, J ) 150 CONTINUE DO 170, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 160, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE ELSE DO 220, J = 1, N TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 190, I = 1, M B( I, J ) = TEMP*B( I, J ) 190 CONTINUE DO 210, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 200, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 200 CONTINUE END IF 210 CONTINUE 220 CONTINUE END IF ELSE * * Form B := alpha*B*A'. * IF( UPPER )THEN DO 260, K = 1, N DO 240, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = ALPHA*A( J, K ) DO 230, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 230 CONTINUE END IF 240 CONTINUE TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) IF( TEMP.NE.ONE )THEN DO 250, I = 1, M B( I, K ) = TEMP*B( I, K ) 250 CONTINUE END IF 260 CONTINUE ELSE DO 300, K = N, 1, -1 DO 280, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = ALPHA*A( J, K ) DO 270, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 270 CONTINUE END IF 280 CONTINUE TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) IF( TEMP.NE.ONE )THEN DO 290, I = 1, M B( I, K ) = TEMP*B( I, K ) 290 CONTINUE END IF 300 CONTINUE END IF END IF END IF * RETURN * * End of DTRMM . * END SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER*4 INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ) * .. * * Purpose * ======= * * DTRMV performs one of the matrix-vector operations * * x := A*x, or x := A'*x, * * where x is an n element vector and A is an n by n unit, or non-unit, * upper or lower triangular matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' x := A*x. * * TRANS = 'T' or 't' x := A'*x. * * TRANS = 'C' or 'c' x := A'*x. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, n ). * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. On exit, X is overwritten with the * tranformed vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER*4 I, INFO, IX, J, JX, KX LOGICAL NOUNIT * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LBAME( UPLO , 'U' ).AND. $ .NOT.LBAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LBAME( TRANS, 'N' ).AND. $ .NOT.LBAME( TRANS, 'T' ).AND. $ .NOT.LBAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LBAME( DIAG , 'U' ).AND. $ .NOT.LBAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOUNIT = LBAME( DIAG, 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( LBAME( TRANS, 'N' ) )THEN * * Form x := A*x. * IF( LBAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = X( J ) DO 10, I = 1, J - 1 X( I ) = X( I ) + TEMP*A( I, J ) 10 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( J, J ) END IF 20 CONTINUE ELSE JX = KX DO 40, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 30, I = 1, J - 1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX + INCX 30 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( J, J ) END IF JX = JX + INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN TEMP = X( J ) DO 50, I = N, J + 1, -1 X( I ) = X( I ) + TEMP*A( I, J ) 50 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( J, J ) END IF 60 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 70, I = N, J + 1, -1 X( IX ) = X( IX ) + TEMP*A( I, J ) IX = IX - INCX 70 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( J, J ) END IF JX = JX - INCX 80 CONTINUE END IF END IF ELSE * * Form x := A'*x. * IF( LBAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 100, J = N, 1, -1 TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 90, I = J - 1, 1, -1 TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE X( J ) = TEMP 100 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 120, J = N, 1, -1 TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 110, I = J - 1, 1, -1 IX = IX - INCX TEMP = TEMP + A( I, J )*X( IX ) 110 CONTINUE X( JX ) = TEMP JX = JX - INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = 1, N TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 130, I = J + 1, N TEMP = TEMP + A( I, J )*X( I ) 130 CONTINUE X( J ) = TEMP 140 CONTINUE ELSE JX = KX DO 160, J = 1, N TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 150, I = J + 1, N IX = IX + INCX TEMP = TEMP + A( I, J )*X( IX ) 150 CONTINUE X( JX ) = TEMP JX = JX + INCX 160 CONTINUE END IF END IF END IF * RETURN * * End of DTRMV . * END SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER*4 M, N, LDA, LDB DOUBLE PRECISION ALPHA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * Purpose * ======= * * DTRSM solves one of the matrix equations * * op( A )*X = alpha*B, or X*op( A ) = alpha*B, * * where alpha is a scalar, X and B are m by n matrices, A is a unit, or * non-unit, upper or lower triangular matrix and op( A ) is one of * * op( A ) = A or op( A ) = A'. * * The matrix X is overwritten on B. * * Parameters * ========== * * SIDE - CHARACTER*1. * On entry, SIDE specifies whether op( A ) appears on the left * or right of X as follows: * * SIDE = 'L' or 'l' op( A )*X = alpha*B. * * SIDE = 'R' or 'r' X*op( A ) = alpha*B. * * Unchanged on exit. * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix A is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANSA - CHARACTER*1. * On entry, TRANSA specifies the form of op( A ) to be used in * the matrix multiplication as follows: * * TRANSA = 'N' or 'n' op( A ) = A. * * TRANSA = 'T' or 't' op( A ) = A'. * * TRANSA = 'C' or 'c' op( A ) = A'. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit triangular * as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of B. M must be at * least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of B. N must be * at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. When alpha is * zero then A is not referenced and B need not be set before * entry. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. * Before entry with UPLO = 'U' or 'u', the leading k by k * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading k by k * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When SIDE = 'L' or 'l' then * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' * then LDA must be at least max( 1, n ). * Unchanged on exit. * * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). * Before entry, the leading m by n part of the array B must * contain the right-hand side matrix B, and on exit is * overwritten by the solution matrix X. * * LDB - INTEGER. * On entry, LDB specifies the first dimension of B as declared * in the calling (sub) program. LDB must be at least * max( 1, m ). * Unchanged on exit. * * * Level 3 Blas routine. * * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LBAME EXTERNAL LBAME * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER*4 I, INFO, J, K, NROWA DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * LSIDE = LBAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LBAME( DIAG , 'N' ) UPPER = LBAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LBAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LBAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LBAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LBAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LBAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LBAME( DIAG , 'U' ) ).AND. $ ( .NOT.LBAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRSM ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * Start the operations. * IF( LSIDE )THEN IF( LBAME( TRANSA, 'N' ) )THEN * * Form B := alpha*inv( A )*B. * IF( UPPER )THEN DO 60, J = 1, N IF( ALPHA.NE.ONE )THEN DO 30, I = 1, M B( I, J ) = ALPHA*B( I, J ) 30 CONTINUE END IF DO 50, K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 40, I = 1, K - 1 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 40 CONTINUE END IF 50 CONTINUE 60 CONTINUE ELSE DO 100, J = 1, N IF( ALPHA.NE.ONE )THEN DO 70, I = 1, M B( I, J ) = ALPHA*B( I, J ) 70 CONTINUE END IF DO 90 K = 1, M IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * Form B := alpha*inv( A' )*B. * IF( UPPER )THEN DO 130, J = 1, N DO 120, I = 1, M TEMP = ALPHA*B( I, J ) DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 120 CONTINUE 130 CONTINUE ELSE DO 160, J = 1, N DO 150, I = M, 1, -1 TEMP = ALPHA*B( I, J ) DO 140, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 140 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF( LBAME( TRANSA, 'N' ) )THEN * * Form B := alpha*B*inv( A ). * IF( UPPER )THEN DO 210, J = 1, N IF( ALPHA.NE.ONE )THEN DO 170, I = 1, M B( I, J ) = ALPHA*B( I, J ) 170 CONTINUE END IF DO 190, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN DO 180, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 180 CONTINUE END IF 190 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 200, I = 1, M B( I, J ) = TEMP*B( I, J ) 200 CONTINUE END IF 210 CONTINUE ELSE DO 260, J = N, 1, -1 IF( ALPHA.NE.ONE )THEN DO 220, I = 1, M B( I, J ) = ALPHA*B( I, J ) 220 CONTINUE END IF DO 240, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN DO 230, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 230 CONTINUE END IF 240 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 250, I = 1, M B( I, J ) = TEMP*B( I, J ) 250 CONTINUE END IF 260 CONTINUE END IF ELSE * * Form B := alpha*B*inv( A' ). * IF( UPPER )THEN DO 310, K = N, 1, -1 IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE END IF DO 290, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 280, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 280 CONTINUE END IF 290 CONTINUE IF( ALPHA.NE.ONE )THEN DO 300, I = 1, M B( I, K ) = ALPHA*B( I, K ) 300 CONTINUE END IF 310 CONTINUE ELSE DO 360, K = 1, N IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 320, I = 1, M B( I, K ) = TEMP*B( I, K ) 320 CONTINUE END IF DO 340, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 330, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 330 CONTINUE END IF 340 CONTINUE IF( ALPHA.NE.ONE )THEN DO 350, I = 1, M B( I, K ) = ALPHA*B( I, K ) 350 CONTINUE END IF 360 CONTINUE END IF END IF END IF * RETURN * * End of DTRSM . * END c ooooooooooooooo SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. Scalar Arguments .. real*8 ALPHA INTEGER*4 INCX, INCY, LDA, M, N * .. Array Arguments .. real*8 A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGER performs the rank 1 operation * * A := alpha*x*y' + A, * * where alpha is a scalar, x is an m element vector, y is an n element * vector and A is an m by n matrix. * * Parameters * ========== * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - real*8. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - real*8 array of dimension at least * ( 1 + ( m - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the m * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - real*8 array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * A - real*8 array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. On exit, A is * overwritten by the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. real*8 ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. real*8 TEMP INTEGER*4 I, INFO, IX, J, JY, KX * .. External Subroutines .. EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( M.LT.0 )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGER ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( INCY.GT.0 )THEN JY = 1 ELSE JY = 1 - ( N - 1 )*INCY END IF IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) DO 10, I = 1, M A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE END IF JY = JY + INCY 20 CONTINUE ELSE IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( M - 1 )*INCX END IF DO 40, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) IX = KX DO 30, I = 1, M A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE END IF JY = JY + INCY 40 CONTINUE END IF * RETURN * * End of DGER . * END SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1992 * * .. Scalar Arguments .. INTEGER*4 INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER*4 IPIV( * ) real*8 A( LDA, * ) * .. * * Purpose * ======= * * DGETF2 computes an LU factorization of a general m-by-n matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 2 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * A (input/output) real*8 array, dimension (LDA,N) * On entry, the m by n matrix to be factored. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER*4 array, dimension (min(M,N)) * The pivot indices; for 1 <= i <= min(M,N), row i of the * matrix was interchanged with row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -k, the k-th argument had an illegal value * > 0: if INFO = k, U(k,k) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations. * * ===================================================================== * * .. Parameters .. real*8 ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER*4 J, JP * .. * .. External Functions .. INTEGER*4 IDAMAX EXTERNAL IDAMAX * .. * .. External Subroutines .. EXTERNAL DGER, DSCAL, DSWAP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETF2', -INFO ) RETURN END IF * * Quick return if possible * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * DO 10 J = 1, MIN( M, N ) * * Find pivot and test for singularity. * JP = J - 1 + IDAMAX( M-J+1, A( J, J ), 1 ) IPIV( J ) = JP IF( A( JP, J ).NE.ZERO ) THEN * * Apply the interchange to columns 1:N. * IF( JP.NE.J ) $ CALL DSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA ) * * Compute elements J+1:M of J-th column. * IF( J.LT.M ) $ CALL DSCAL( M-J, ONE / A( J, J ), A( J+1, J ), 1 ) * ELSE IF( INFO.EQ.0 ) THEN * INFO = J END IF * IF( J.LT.MIN( M, N ) ) THEN * * Update trailing submatrix. * CALL DGER( M-J, N-J, -ONE, A( J+1, J ), 1, A( J, J+1 ), LDA, $ A( J+1, J+1 ), LDA ) END IF 10 CONTINUE RETURN * * End of DGETF2 * END SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * March 31, 1993 * * .. Scalar Arguments .. INTEGER*4 INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER*4 IPIV( * ) real*8 A( LDA, * ) * .. * * Purpose * ======= * * DGETRF computes an LU factorization of a general M-by-N matrix A * using partial pivoting with row interchanges. * * The factorization has the form * A = P * L * U * where P is a permutation matrix, L is lower triangular with unit * diagonal elements (lower trapezoidal if m > n), and U is upper * triangular (upper trapezoidal if m < n). * * This is the right-looking Level 3 BLAS version of the algorithm. * * Arguments * ========= * * M (input) INT * The number of rows of the matrix A. M >= 0. * * N (input) INT * The number of columns of the matrix A. N >= 0. * * A (input/output) real*8 array, dimension (LDA,N) * On entry, the M-by-N matrix to be factored. * On exit, the factors L and U from the factorization * A = P*L*U; the unit diagonal elements of L are not stored. * * LDA (input) INT * The leading dimension of the array A. LDA >= max(1,M). * * IPIV (output) INTEGER*4 array, dimension (min(M,N)) * The pivot indices; for 1 <= i <= min(M,N), row i of the * matrix was interchanged with row IPIV(i). * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, U(i,i) is exactly zero. The factorization * has been completed, but the factor U is exactly * singular, and division by zero will occur if it is used * to solve a system of equations. * * ===================================================================== * * .. Parameters .. real*8 ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER*4 I, IINFO, J, JB, NB * .. * .. External Subroutines .. EXTERNAL DGEMM, DGETF2, DLASWP, DTRSM, XERBLA * .. * .. External Functions .. INTEGER*4 ILAENV EXTERNAL ILAENV * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DGETRF', -INFO ) RETURN END IF * * Quick return if possible * IF( M.EQ.0 .OR. N.EQ.0 ) $ RETURN * * Determine the block size for this environment. * NB = ILAENV( 1, 'DGETRF', ' ', M, N, -1, -1 ) IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN * * Use unblocked code. * CALL DGETF2( M, N, A, LDA, IPIV, INFO ) ELSE * * Use blocked code. * DO 20 J = 1, MIN( M, N ), NB JB = MIN( MIN( M, N )-J+1, NB ) * * Factor diagonal and subdiagonal blocks and test for exact * singularity. * CALL DGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO ) * * Adjust INFO and the pivot indices. * IF( INFO.EQ.0 .AND. IINFO.GT.0 ) $ INFO = IINFO + J - 1 DO 10 I = J, MIN( M, J+JB-1 ) IPIV( I ) = J - 1 + IPIV( I ) 10 CONTINUE * * Apply interchanges to columns 1:J-1. * CALL DLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 ) * IF( J+JB.LE.N ) THEN * * Apply interchanges to columns J+JB:N. * CALL DLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, $ IPIV, 1 ) * * Compute block row of U. * CALL DTRSM( 'L', 'L', 'N', 'U', JB, $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ), $ LDA ) IF( J+JB.LE.M ) THEN * * Update trailing submatrix. * CALL DGEMM( 'N', 'N', M-J-JB+1, $ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA, $ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ), $ LDA ) END IF END IF 20 CONTINUE END IF RETURN * * End of DGETRF * END SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1999 * * .. Scalar Arguments .. INTEGER*4 INCX, K1, K2, LDA, N * .. * .. Array Arguments .. INTEGER*4 IPIV( * ) real*8 A( LDA, * ) * .. * * Purpose * ======= * * DLASWP performs a series of row interchanges on the matrix A. * One row interchange is initiated for each of rows K1 through K2 of A. * * Arguments * ========= * * N (input) INTEGER * The number of columns of the matrix A. * * A (input/output) real*8 array, dimension (LDA,N) * On entry, the matrix of column dimension N to which the row * interchanges will be applied. * On exit, the permuted matrix. * * LDA (input) INTEGER * The leading dimension of the array A. * * K1 (input) INTEGER * The first element of IPIV for which a row interchange will * be done. * * K2 (input) INTEGER * The last element of IPIV for which a row interchange will * be done. * * IPIV (input) INTEGER*4 array, dimension (M*abs(INCX)) * The vector of pivot indices. Only the elements in positions * K1 through K2 of IPIV are accessed. * IPIV(K) = L implies rows K and L are to be interchanged. * * INCX (input) INTEGER * The increment between successive values of IPIV. If IPIV * is negative, the pivots are applied in reverse order. * * Further Details * =============== * * Modified by * R. C. Whaley, Computer Science Dept., Univ. of Tenn., Knoxville, USA * * ===================================================================== * * .. Local Scalars .. INTEGER*4 I, I1, I2, INC, IP, IX, IX0, J, K, N32 real*8 TEMP * .. * .. Executable Statements .. * * Interchange row I with row IPIV(I) for each of rows K1 through K2. * IF( INCX.GT.0 ) THEN IX0 = K1 I1 = K1 I2 = K2 INC = 1 ELSE IF( INCX.LT.0 ) THEN IX0 = 1 + ( 1-K2 )*INCX I1 = K2 I2 = K1 INC = -1 ELSE RETURN END IF * N32 = ( N / 32 )*32 IF( N32.NE.0 ) THEN DO 30 J = 1, N32, 32 IX = IX0 DO 20 I = I1, I2, INC IP = IPIV( IX ) IF( IP.NE.I ) THEN DO 10 K = J, J + 31 TEMP = A( I, K ) A( I, K ) = A( IP, K ) A( IP, K ) = TEMP 10 CONTINUE END IF IX = IX + INCX 20 CONTINUE 30 CONTINUE END IF IF( N32.NE.N ) THEN N32 = N32 + 1 IX = IX0 DO 50 I = I1, I2, INC IP = IPIV( IX ) IF( IP.NE.I ) THEN DO 40 K = N32, N TEMP = A( I, K ) A( I, K ) = A( IP, K ) A( IP, K ) = TEMP 40 CONTINUE END IF IX = IX + INCX 50 CONTINUE END IF * RETURN * * End of DLASWP * 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 modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c real*8 da,dx(*) integer*4 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 integer*4 function idamax(n,dx,incx) c c finds the index of element having max. absolute value. 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 real*8 dx(*),dmax integer*4 i,incx,ix,n c idamax = 0 if( n.lt.1 .or. incx.le.0 ) 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 INTEGER*4 FUNCTION IEEECK( ISPEC, ZERO, ONE ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1998 * * .. Scalar Arguments .. INTEGER*4 ISPEC REAL*8 ONE, ZERO * .. * * Purpose * ======= * * IEEECK is called from the ILAENV to verify that Infinity and * possibly NaN arithmetic is safe (i.e. will not trap). * * Arguments * ========= * * ISPEC (input) INTEGER * Specifies whether to test just for inifinity arithmetic * or whether to test for infinity and NaN arithmetic. * = 0: Verify infinity arithmetic only. * = 1: Verify infinity and NaN arithmetic. * * ZERO (input) REAL * Must contain the value 0.0 * This is passed to prevent the compiler from optimizing * away this code. * * ONE (input) REAL * Must contain the value 1.0 * This is passed to prevent the compiler from optimizing * away this code. * * RETURN VALUE: INTEGER * = 0: Arithmetic failed to produce the correct answers * = 1: Arithmetic produced the correct answers * * .. Local Scalars .. REAL*8 NAN1, NAN2, NAN3, NAN4, NAN5, NAN6, NEGINF, $ NEGZRO, NEWZRO, POSINF * .. * .. Executable Statements .. IEEECK = 1 * POSINF = ONE / ZERO IF( POSINF.LE.ONE ) THEN IEEECK = 0 RETURN END IF * NEGINF = -ONE / ZERO IF( NEGINF.GE.ZERO ) THEN IEEECK = 0 RETURN END IF * NEGZRO = ONE / ( NEGINF+ONE ) IF( NEGZRO.NE.ZERO ) THEN IEEECK = 0 RETURN END IF * NEGINF = ONE / NEGZRO IF( NEGINF.GE.ZERO ) THEN IEEECK = 0 RETURN END IF * NEWZRO = NEGZRO + ZERO IF( NEWZRO.NE.ZERO ) THEN IEEECK = 0 RETURN END IF * POSINF = ONE / NEWZRO IF( POSINF.LE.ONE ) THEN IEEECK = 0 RETURN END IF * NEGINF = NEGINF*POSINF IF( NEGINF.GE.ZERO ) THEN IEEECK = 0 RETURN END IF * POSINF = POSINF*POSINF IF( POSINF.LE.ONE ) THEN IEEECK = 0 RETURN END IF * * * * * Return if we were only asked to check infinity arithmetic * IF( ISPEC.EQ.0 ) $ RETURN * NAN1 = POSINF + NEGINF * NAN2 = POSINF / NEGINF * NAN3 = POSINF / POSINF * NAN4 = POSINF*ZERO * NAN5 = NEGINF*NEGZRO * NAN6 = NAN5*0.0d0 * IF( NAN1.EQ.NAN1 ) THEN IEEECK = 0 RETURN END IF * IF( NAN2.EQ.NAN2 ) THEN IEEECK = 0 RETURN END IF * IF( NAN3.EQ.NAN3 ) THEN IEEECK = 0 RETURN END IF * IF( NAN4.EQ.NAN4 ) THEN IEEECK = 0 RETURN END IF * IF( NAN5.EQ.NAN5 ) THEN IEEECK = 0 RETURN END IF * IF( NAN6.EQ.NAN6 ) THEN IEEECK = 0 RETURN END IF * RETURN END INTEGER*4 FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, $ N4 ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * June 30, 1999 * * .. Scalar Arguments .. CHARACTER*( * ) NAME, OPTS INTEGER*4 ISPEC, N1, N2, N3, N4 * .. * * Purpose * ======= * * ILAENV is called from the LAPACK routines to choose problem-dependent * parameters for the local environment. See ISPEC for a description of * the parameters. * * This version provides a set of parameters which should give good, * but not optimal, performance on many of the currently available * computers. Users are encouraged to modify this subroutine to set * the tuning parameters for their particular machine using the option * and problem size information in the arguments. * * This routine will not function correctly if it is converted to all * lower case. Converting it to all upper case is allowed. * * Arguments * ========= * * ISPEC (input) INTEGER * Specifies the parameter to be returned as the value of * ILAENV. * = 1: the optimal blocksize; if this value is 1, an unblocked * algorithm will give the best performance. * = 2: the minimum block size for which the block routine * should be used; if the usable block size is less than * this value, an unblocked routine should be used. * = 3: the crossover point (in a block routine, for N less * than this value, an unblocked routine should be used) * = 4: the number of shifts, used in the nonsymmetric * eigenvalue routines * = 5: the minimum column dimension for blocking to be used; * rectangular blocks must have dimension at least k by m, * where k is given by ILAENV(2,...) and m by ILAENV(5,...) * = 6: the crossover point for the SVD (when reducing an m by n * matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds * this value, a QR factorization is used first to reduce * the matrix to a triangular form.) * = 7: the number of processors * = 8: the crossover point for the multishift QR and QZ methods * for nonsymmetric eigenvalue problems. * = 9: maximum size of the subproblems at the bottom of the * computation tree in the divide-and-conquer algorithm * (used by xGELSD and xGESDD) * =10: ieee NaN arithmetic can be trusted not to trap * =11: infinity arithmetic can be trusted not to trap * * NAME (input) CHARACTER*(*) * The name of the calling subroutine, in either upper case or * lower case. * * OPTS (input) CHARACTER*(*) * The character options to the subroutine NAME, concatenated * into a single character string. For example, UPLO = 'U', * TRANS = 'T', and DIAG = 'N' for a triangular routine would * be specified as OPTS = 'UTN'. * * N1 (input) INTEGER * N2 (input) INTEGER * N3 (input) INTEGER * N4 (input) INTEGER * Problem dimensions for the subroutine NAME; these may not all * be required. * * (ILAENV) (output) INTEGER * >= 0: the value of the parameter specified by ISPEC * < 0: if ILAENV = -k, the k-th argument had an illegal value. * * Further Details * =============== * * The following conventions have been used when calling ILAENV from the * LAPACK routines: * 1) OPTS is a concatenation of all of the character options to * subroutine NAME, in the same order that they appear in the * argument list for NAME, even if they are not used in determining * the value of the parameter specified by ISPEC. * 2) The problem dimensions N1, N2, N3, N4 are specified in the order * that they appear in the argument list for NAME. N1 is used * first, N2 second, and so on, and unused problem dimensions are * passed a value of -1. * 3) The parameter value returned by ILAENV is checked for validity in * the calling subroutine. For example, ILAENV is used to retrieve * the optimal blocksize for STRTRI as follows: * * NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) * IF( NB.LE.1 ) NB = MAX( 1, N ) * * ===================================================================== * * .. Local Scalars .. LOGICAL CNAME, SNAME CHARACTER*1 C1 CHARACTER*2 C2, C4 CHARACTER*3 C3 CHARACTER*6 SUBNAM INTEGER*4 I, IC, IZ, NB, NBMIN, NX * .. * .. Intrinsic Functions .. INTRINSIC CHAR, ICHAR, INT, MIN, REAL * .. * .. External Functions .. INTEGER*4 IEEECK EXTERNAL IEEECK * .. * .. Executable Statements .. * GO TO ( 100, 100, 100, 400, 500, 600, 700, 800, 900, 1000, $ 1100 ) ISPEC * * Invalid value for ISPEC * ILAENV = -1 RETURN * 100 CONTINUE * * Convert NAME to upper case if the first character is lower case. * ILAENV = 1 SUBNAM = NAME IC = ICHAR( SUBNAM( 1:1 ) ) IZ = ICHAR( 'Z' ) IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN * * ASCII character set * IF( IC.GE.97 .AND. IC.LE.122 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 10 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.97 .AND. IC.LE.122 ) $ SUBNAM( I:I ) = CHAR( IC-32 ) 10 CONTINUE END IF * ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN * * EBCDIC character set * IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. $ ( IC.GE.162 .AND. IC.LE.169 ) ) THEN SUBNAM( 1:1 ) = CHAR( IC+64 ) DO 20 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR. $ ( IC.GE.145 .AND. IC.LE.153 ) .OR. $ ( IC.GE.162 .AND. IC.LE.169 ) ) $ SUBNAM( I:I ) = CHAR( IC+64 ) 20 CONTINUE END IF * ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN * * Prime machines: ASCII+128 * IF( IC.GE.225 .AND. IC.LE.250 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 30 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.225 .AND. IC.LE.250 ) $ SUBNAM( I:I ) = CHAR( IC-32 ) 30 CONTINUE END IF END IF * C1 = SUBNAM( 1:1 ) SNAME = C1.EQ.'S' .OR. C1.EQ.'D' CNAME = C1.EQ.'C' .OR. C1.EQ.'Z' IF( .NOT.( CNAME .OR. SNAME ) ) $ RETURN C2 = SUBNAM( 2:3 ) C3 = SUBNAM( 4:6 ) C4 = C3( 2:3 ) * GO TO ( 110, 200, 300 ) ISPEC * 110 CONTINUE * * ISPEC = 1: block size * * In these examples, separate code is provided for setting NB for * real and complex. We assume that NB will take the same value in * single or real*8. * NB = 1 * IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'PO' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NB = 32 ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRF' ) THEN NB = 64 ELSE IF( C3.EQ.'TRD' ) THEN NB = 32 ELSE IF( C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( C2.EQ.'GB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'PB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'TR' ) THEN IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'LA' ) THEN IF( C3.EQ.'UUM' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN IF( C3.EQ.'EBZ' ) THEN NB = 1 END IF END IF ILAENV = NB RETURN * 200 CONTINUE * * ISPEC = 2: minimum block size * NBMIN = 2 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NBMIN = 8 ELSE NBMIN = 8 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF END IF ILAENV = NBMIN RETURN * 300 CONTINUE * * ISPEC = 3: crossover point * NX = 0 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR. $ C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( SNAME .AND. C3.EQ.'TRD' ) THEN NX = 32 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NX = 32 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NX = 128 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR. $ C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR. $ C4.EQ.'BR' ) THEN NX = 128 END IF END IF END IF ILAENV = NX RETURN * 400 CONTINUE * * ISPEC = 4: number of shifts (used by xHSEQR) * ILAENV = 6 RETURN * 500 CONTINUE * * ISPEC = 5: minimum column dimension (not used) * ILAENV = 2 RETURN * 600 CONTINUE * * ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) * ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 ) RETURN * 700 CONTINUE * * ISPEC = 7: number of processors (not used) * ILAENV = 1 RETURN * 800 CONTINUE * * ISPEC = 8: crossover point for multishift (used by xHSEQR) * ILAENV = 50 RETURN * 900 CONTINUE * * ISPEC = 9: maximum size of the subproblems at the bottom of the * computation tree in the divide-and-conquer algorithm * (used by xGELSD and xGESDD) * ILAENV = 25 RETURN * 1000 CONTINUE * * ISPEC = 10: ieee NaN arithmetic can be trusted not to trap * C ILAENV = 0 ILAENV = 1 IF( ILAENV.EQ.1 ) THEN ILAENV = IEEECK( 0, 0.0D0, 1.0D0 ) END IF RETURN * 1100 CONTINUE * * ISPEC = 11: infinity arithmetic can be trusted not to trap * C ILAENV = 0 ILAENV = 1 IF( ILAENV.EQ.1 ) THEN ILAENV = IEEECK( 1, 0.0D0, 1.0D0 ) END IF RETURN * * End of ILAENV * END c Updated 10/24/2001. c ccccccccccccccccccccccccc Program 4.2 cccccccccccccccccccccccccc c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c Please Note: c c c c (1) This computer program is part of the book, "An Introduction to c c Computational Physics," written by Tao Pang and published and c c copyrighted by Cambridge University Press in 1997. c c c c (2) No warranties, express or implied, are made for this program. c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c SUBROUTINE DTRM(A,M,N,D,INDX,C) C C Subroutine for evaluating the determinant of a matrix using C the partial-pivoting Gaussian elimination scheme. C implicit none integer*4 M,N,INDX(M),I,MSGN,J real*8 A(M,M),D,C(M) C CALL ELGS(A,M,N,INDX,C) D = 1.0d0 DO 100 I = 1, N 100 D = D*A(INDX(I),I) MSGN = 1 DO 200 I = 1, N DO 150 WHILE (I.NE.INDX(I)) MSGN = -MSGN J = INDX(I) INDX(I) = INDX(J) INDX(J) = J 150 END DO 200 CONTINUE D = dble(MSGN)*D RETURN END C SUBROUTINE ELGS(A,M,N,INDX,C) C C Subroutine to perform the partial-pivoting Gaussian elimination. C A(N,N) is the original matrix in the input and transformed C matrix plus the pivoting element ratios below the diagonal in C the output. INDX(N) records the pivoting order. C implicit none integer*4 M,N,INDX(M),I,J,K,ITMP real*8 A(M,M),C(M),C1,PI1,PI,PJ C DO 50 I = 1, N 50 INDX(I) = I C C Find the rescaling factors, one from each row C DO 100 I = 1, N C1= 0.0d0 DO 90 J = 1, N 90 C1 = MAX(C1,DABS(A(I,J))) 100 C(I) = C1 C C Search the pivoting (largest) element from each column C DO 200 J = 1, N-1 PI1 = 0.0d0 DO 150 I = J, N PI = ABS(A(INDX(I),J))/C(INDX(I)) IF (PI.GT.PI1) THEN PI1 = PI K = I ELSE ENDIF 150 CONTINUE C C Interchange the rows via INDX(N) to record pivoting order C ITMP = INDX(J) INDX(J) = INDX(K) INDX(K) = ITMP DO 170 I = J+1, N PJ = A(INDX(I),J)/A(INDX(J),J) C C Record pivoting ratios below the diagonal C A(INDX(I),J) = PJ C C Modify other elements accordingly C DO 160 K = J+1, N A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K) 160 CONTINUE 170 CONTINUE 200 CONTINUE C RETURN END