program genchar c c generates M charges around a molecule of N atoms c c CHAR.OPT ... file with options c NAMES.LST ... list of MD geometries c c NMC - number of MD confugurations c LPC - (logical) - bondary conditions c LAV - (logical) - construct average charge geometries c LAM - (logical) - construct average solvent molecular geometries c LAT - use MD (Tinker) charges for average solvent c LIN - interactive run c LSP - get potential from scratch POT.SCR c LSX - get potential from text file POT.TXT c LQ0 - minimize RMS from Tinker charges c LWR - larger output c DPCX,DPCY,DPCZ - lengths of box c N1,N2 - the potential is calculated at atoms N1...N2,NAT=N2-N1+1 c NC - number of charges, NC.GE.NAT+2 c fnpar - filename with MD force field c NAM0 - maximum number of atoms in one solvent molecule c maxit - maximum number of iteration for the group division c tol - maximum allowed change in the moment sum c M0 max number of atoms in one cube/MD configuration c implicit real*8 (a-h,o-z) implicit integer*4 (i-n) parameter (nat0=400,nc0=1000,M0=1000,nc00=1000, 1nt0=4,NAM0=10,natall=500000) character*8 nte character*1 s1 character*80 fn,s80,fnpar logical LPC,LWR,LAV,LAM,LAT,LIN,LQ0,LSP,LSX dimension ct(nc0),x(nat0),y(nat0),z(nat0),cti(M0), 1izt(nc0),izq(nat0),c(M0),b(M0),a(M0,M0),ipiv(M0), 2work(3*M0*M0),q(nat0),xi(nc00),yi(nc00),zi(nc00), 3qi(nc00),di(nc00),a0(M0,M0),itk(nt0),ntk(nt0),xit(nc00), 4yit(nc00),zit(nc00),alistx(nc00,NAM0),iti(nc00), 5alisty(nc00,NAM0),alistz(nc00,NAM0),listi(nc00,NAM0) real*4 rijx(natall),rijy(natall),rijz(natall), 1tkx(nc0),tky(nc0),tkz(nc0) common/bigc/rijx,rijy,rijz dimension xo(NAM0),yo(NAM0),zo(NAM0),xn(NAM0),yn(NAM0),zn(NAM0), 1xh(NAM0),yh(NAM0),zh(NAM0),iatom(natall), 2tlistx(natall,NAM0),tlisty(natall,NAM0),tlistz(natall,NAM0), 3itlist(natall,NAM0) common/mv2/xo,yo,zo,xh,yh,zh,tlistx,tlisty,tlistz,xn,yn,zn, 1NAM,NMC,ii1,iatom,itlist common/inter/LIN c tol=0.0d0 maxit=2000 LWR=.false. LAV=.false. LAM=.false. LAT=.false. LPC=.false. LSP=.false. LSX=.false. LQ0=.false. LIN=.true. NMC=0 NC=0 N1=1 n2=0 NAK=2 NAM=3 itk(1)=8 itk(2)=1 open(3,file='CHAR.OUT') call iopar(LAV,LAM,LAT,LPC,LSP,LQ0,LWR,LIN,NMC,NC, 1DPCX,DPCY,DPCZ,NAK,NAM,itk,ntk,ntm,N1,N2,fnpar, 2maxit,tol,toi,LSX) NAT=N2-N1+1 if(LIN)write(6,*)NC,NAT write(3,*)NC,NAT if(NC.lt.NAT)call report('NC < NAT !') if(NC.gt.nc00)call report('NC > nc00 !') if(NAT.gt.nat0)call report('nat0 < NAT !') if(NAM.gt.NAM0)call report('NAM < NAM0 !') if(NC.ne.NTM*NAM.and.LAM)call report('NC <> NAM * NTM !') M=NC+NAT+1 if(M.gt.M0)call report('M > M0!') call co(NMC,NC,LPC,DPCX,DPCY,DPCZ,N1,N2,fnpar,NAT,M) do 43 i=1,M 43 c(i)=0.0d0 dx2=DPCX/2.0d0 dy2=DPCY/2.0d0 dz2=DPCZ/2.0d0 c c assign charges and atomic numbers to types do 7 i=1,nc0 izt(i)=0 7 ct(i)=0.0d0 nt=0 nz=0 open(2,file=fnpar,status='OLD') 8 read(2,2828,end=9,err=9)s80 2828 format(a80) if(s80(1:4).eq.'atom')then open(55,file='scr') write(55,2828)s80(5:12)//s80(55:70) rewind(55) read(55,*)it,izta if(it.lt.nc0)then izt(it)=izta nz=nz+1 c write(6,6007)nz,it,izt(it) c6007 format(3i5) else call report('Too many atom types') endif close(55) endif if(s80(1:6).eq.'charge')then open(55,file='scr') write(55,2828)s80(7:80) rewind(55) read(55,*)it,qt if(it.lt.nc0)then ct(it)=qt nt=nt+1 c write(6,6006)nt,it,qt c6006 format(2i5,f8.4) else call report('Too many charge types') endif close(55) endif goto 8 9 close(2) write(3,*)nt,' charge types found in '//fnpar(1:40) if(LIN)write(6,*)nt,' charge types found in '//fnpar(1:40) write(3,*)nz,' atom types found in '//fnpar(1:40) if(LIN)write(6,*)nz,' atom types found in '//fnpar(1:40) open(2,file='NAMES.LST',status='OLD') do 3 i=1,NMC read(2,300)fn 300 format(a80) c write(6,*)fn(1:60) open(21,file=fn,status='OLD') call ratm(N1,N2,x,y,z,q,ct,izq,izt,nc0,LSX) rewind 21 read(21,*)natmd do 6 ia=1,natmd read(21,2121)s1,xa,ya,za,it 2121 format(6x,2x,a1,2x,3f12.6,i6) if(ia.lt.N1.or.ia.gt.N2)then if(it.eq.0.and.LSX)then if(s1.eq.'O')it=649 if(s1.eq.'H')it=650 c N,C,glycin types in amber.prm: if(s1.eq.'C')it=2 if(s1.eq.'N')it=1 izt(1)=7 izt(2)=6 izt(649)=8 izt(650)=1 endif if(it.lt.1.or.it.gt.nc0)call report('Invalid it.') qc=ct(it) do 10 ja=1,NAT xj=x(ja) yj=y(ja) zj=z(ja) if(LPC)then if(xa-xj.gt.dx2)xa=xa-DPCX if(xj-xa.gt.dx2)xa=xa+DPCX if(ya-yj.gt.dy2)ya=ya-DPCY if(yj-ya.gt.dy2)ya=ya+DPCY if(za-zj.gt.dz2)za=za-DPCZ if(zj-za.gt.dz2)za=za+DPCZ endif rc=dsqrt((xj-xa)**2+(yj-ya)**2+(zj-za)**2) 10 c(nc+ja)=c(nc+ja)+qc/rc endif 6 continue close(21) 3 continue write(3,*)NMC,' configurations read in' if(LIN)write(6,*)NMC,' configurations read in' close(2) if(LSX)then icr=0 icl=0 open(8,file='POT.TXT',status='old') 2011 read(8,808,end=333,err=333)s1 808 format(a1) icl=icl+1 if(s1.ne.' ')then backspace(8) read(8,*)ja,idum,idum,c(nc+ja) icr=icr+1 endif goto 2011 333 close(8) write(3,*)'POT.TXT read,',icl,' lines,',icr,' charges.' if(LIN)write(6,*)'POT.TXT read,',icl,' lines,',icr,' charges.' endif if(LSP)then write(3,*)'POT.SCR read' if(LIN)write(6,*)'POT.SCR read' else if(LIN)write(6,*)'Average potentials' write(3,*)'Average potentials' endif open(8,file='POT.SCR',form='unformatted') do 11 ja=1,NAT if(LSP)then read(8)c(nc+ja) else c(nc+ja)=c(nc+ja)/dble(NMC) write(8)c(nc+ja) endif if(LIN)write(6,611)ja,izq(ja),c(nc+ja) 11 write(3,611)ja,izq(ja),c(nc+ja) 611 format(2i6,f10.4) close(8) c c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(LAV)then c c find "average" positions for solvent atoms c NAK ... Number of atoms kinds (e.g. 2 for h2o) c itk ... list of atom numbers c ntk ... number of each kind c c Loop for each kind of atoms ii=0 do 1001 ik=1,NAK write(3,*)'Selecting kind ',itk(ik),',',ntk(ik),' atoms' if(LIN)write(6,*)'Selecting kind ',itk(ik),',',ntk(ik),' atoms' ito=0 open(2,file='NAMES.LST',status='OLD') do 1003 i=1,NMC read(2,300)fn open(21,file=fn,status='OLD') call ratm(N1,N2,x,y,z,q,ct,izq,izt,nc0,LSX) rewind 21 read(21,*)natmd ic=0 do 161 ia=1,natmd read(21,2121)s1,xa,ya,za,it if(it.eq.0.and.LSX)then if(s1.eq.'O')it=649 if(s1.eq.'H')it=650 c N,C,glycin types in amber.prm: if(s1.eq.'C')it=2 if(s1.eq.'N')it=1 izt(1)=7 izt(2)=6 izt(649)=8 izt(650)=1 endif if(ia.lt.N1.or.ia.gt.N2.and.it.eq.itk(ik))then qc=ct(it) ic=ic+1 xit(ic)=xa yit(ic)=ya zit(ic)=za di(ic)=dmin(NAT,x,y,z,LPC,xa,ya,za,dx2,dy2,dz2,DPCX, 1DPCY,DPCZ) qi(ic)=qc endif 161 continue close(21) if(ic.lt.ntk(ik))call report('Not enough charges was found') c c order charges do 281 ir=1,ic do 281 jr=ir+1,ic 281 if(di(jr).lt.di(ir))call sw(ir,jr,xit,yit,zit,qi,di) c c take ntk(ik) of them do 1004 ir=1,ntk(ik) ito=ito+1 if(ito.gt.natall)call report('ito overflow') rijx(ito)=real(xit(ir)) rijy(ito)=real(yit(ir)) 1004 rijz(ito)=real(zit(ir)) 1003 continue c c now make NMC groups, each with ntk(ik) atoms open(45,file='tem1.x') write(45,*) write(45,*)ito do 10051 i=1,ito write(45,41)1,rijx(i),rijy(i),rijz(i) 10051 continue close(45) write(3,*)'tem1',ik if(LIN)write(6,*)'tem1',ik if(LIN)read(5,*)idumm call makegr(NMC,ntk(ik),tkx,tky,tkz,iatom,maxit,tol) c c record the average charge positions open(45,file='tem2.x') write(45,*) write(45,*)ntk(ik) do 1005 i=1,ntk(ik) ii=ii+1 xi(ii)=dble(tkx(i)) yi(ii)=dble(tky(i)) zi(ii)=dble(tkz(i)) write(45,41)1,xi(ii),yi(ii),zi(ii) 41 format(i5,3f12.6) 1005 continue close(45) write(3,*)'tem2' if(LIN)write(6,*)'tem2' if(LIN)read(5,*)idumm close(2) 1001 continue else c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c ***************************************************************** if(LAM)then c c find "average" positions for solvent molecules c NAM ... Number of atoms in the solvent molecule c ntm ... number of molecules of solvent to be found c write(3,*)'Selecting molecules each with ',NAM,' atoms' if(LIN)write(6,*)'Selecting molecules each with ',NAM,' atoms' ito=0 open(2,file='NAMES.LST',status='OLD') do 10031 i=1,NMC read(2,300)fn open(21,file=fn,status='OLD') c c read solute: call ratm(N1,N2,x,y,z,q,ct,izq,izt,nc0,LSX) rewind 21 c c read again to load molecules imol=0 read(21,*)natmd ia=0 162 ia=ia+1 if(ia.lt.N1.or.ia.gt.N2)then ctx=0.0d0 cty=0.0d0 ctz=0.0d0 imol=imol+1 c c read first atom extra so that periodic conditions can be fixed read(21,2121)s1,xj,yj,zj,it if(it.eq.0.and.LSX)then if(s1.eq.'O')it=649 if(s1.eq.'H')it=650 c N,C,glycin types in amber.prm: if(s1.eq.'C')it=2 if(s1.eq.'N')it=1 izt(1)=7 izt(2)=6 izt(649)=8 izt(650)=1 endif alistx(imol,1)=xj alisty(imol,1)=yj alistz(imol,1)=zj listi(imol,1)=it ctx=ctx+xj cty=cty+yj ctz=ctz+zj do 163 iia=2,NAM read(21,2121)s1,xa,ya,za,it if(it.eq.0.and.LSX)then if(s1.eq.'O')it=649 if(s1.eq.'H')it=650 c N,C,glycin types in amber.prm: if(s1.eq.'C')it=2 if(s1.eq.'N')it=1 izt(1)=7 izt(2)=6 izt(649)=8 izt(650)=1 endif if(LPC)then if(xa-xj.gt.dx2)xa=xa-DPCX if(xj-xa.gt.dx2)xa=xa+DPCX if(ya-yj.gt.dy2)ya=ya-DPCY if(yj-ya.gt.dy2)ya=ya+DPCY if(za-zj.gt.dz2)za=za-DPCZ if(zj-za.gt.dz2)za=za+DPCZ endif alistx(imol,iia)=xa alisty(imol,iia)=ya alistz(imol,iia)=za listi(imol,iia)=it ctx=ctx+xa cty=cty+ya 163 ctz=ctz+za ctx=ctx/dble(NAM) cty=cty/dble(NAM) ctz=ctz/dble(NAM) xit(imol)=ctx yit(imol)=cty zit(imol)=ctz ia=ia+NAM-1 c c determine distance of the molecular center to the solute: di(imol)=dmin(NAT,x,y,z,LPC,ctx,cty,ctz,dx2,dy2,dz2,DPCX, 1 DPCY,DPCZ) else read(21,*) endif if(ia.lt.natmd)goto 162 close(21) if(imol.lt.ntm)call report('Not enough solvent molecules found') c c order solvent do 2811 ir=1,imol do 2811 jr=ir+1,imol 2811 if(di(jr).lt.di(ir))call swm(ir,jr,xit,yit,zit,qi,di, 1alistx,alisty,alistz,NAM0,NAM,listi,nc00) c c take first ntm of them do 10041 ir=1,ntm ito=ito+1 if(ito.gt.natall)call report('ito overflow') do 10042 ia=1,NAM itlist(ito,ia)=listi(ir,ia) tlistx(ito,ia)=alistx(ir,ia) tlisty(ito,ia)=alisty(ir,ia) 10042 tlistz(ito,ia)=alistz(ir,ia) rijx(ito)=real(xit(ir)) rijy(ito)=real(yit(ir)) 10041 rijz(ito)=real(zit(ir)) c c record geometry of this configuration open(8,file='scr') write(8,8000)i,'MC.x' 8000 format(i6,a4) rewind 8 read(8,8001)nte 8001 format(a10) close(8) if(LWR)then do 101 ist=1,len(nte) 101 if(nte(ist:ist).ne.' ')goto 102 102 open(8,file=nte(ist:len(nte))) write(8,8001)nte write(8,*)ntm*NAM+N2-N1+1 do 20041 ia=1,N2-N1+1 20041 write(8,41)izq(ia),x(ia),y(ia),z(ia) do 20042 ir=1,ntm do 20042 ia=1,NAM it=1 if(listi(ir,ia).eq.649)it=8 20042 write(8,41)it,alistx(ir,ia),alisty(ir,ia),alistz(ir,ia) close(8) endif 10031 continue close(2) open(45,file='temallm.x') write(45,*) write(45,*)ito*(1+NAM) do 11051 i=1,ito write(45,41)9,rijx(i),rijy(i),rijz(i) do 11051 j=1,NAM it=1 if(itlist(i,j).eq.649)it=8 11051 write(45,41)it,tlistx(i,j),tlisty(i,j),tlistz(i,j) close(45) c c now make NMC groups, each with ntm molecules call makegr(NMC,ntm,tkx,tky,tkz,iatom,maxit,tol) c c record the average positions open(45,file='avem.x') write(45,*) write(45,*)ntm do 1105 i=1,ntm 1105 write(45,41)9,tkx(i),tky(i),tkz(i) close(45) write(3,*)'avem.x' if(LIN)write(6,*)'avem.x' c open(45,file='whole.x') open(46,file='whole.scr',form='unformatted') write(45,*) write(45,*)ntm*NAM write(46)ntm*NAM c for each molecule, get the best orientation do 1106 im=1,ntm c c get the first molecule in the group ii1=NMC*(im-1)+1 iig=iatom(ii1) tx=0.0d0 ty=0.0d0 tz=0.0d0 do 1107 ia=1,NAM xo(ia)=tlistx(iig,ia) yo(ia)=tlisty(iig,ia) zo(ia)=tlistz(iig,ia) tx=tx+xo(ia) ty=ty+yo(ia) 1107 tz=tz+zo(ia) c c subtract center do 1108 ia=1,NAM xo(ia)=xo(ia)-tx/dble(NAM) yo(ia)=yo(ia)-ty/dble(NAM) 1108 zo(ia)=zo(ia)-tz/dble(NAM) c c optimize orientation call mopti c c shift to center and write to file do 1109 ia=1,NAM xo(ia)=xn(ia)+dble(tkx(im)) yo(ia)=yn(ia)+dble(tky(im)) zo(ia)=zn(ia)+dble(tkz(im)) it=1 if(itlist(iig,ia).eq.649)it=8 write(46)itlist(iig,ia),xo(ia),yo(ia),zo(ia) 1109 write(45,41)it,xo(ia),yo(ia),zo(ia) 1106 continue close(45) close(46) write(3,*)'whole.x' if(LIN)write(6,*)'whole.x' c c calculate average potential with tinker charges do 1110 ja=1,NAT 1110 cti(nc+ja)=0.0d0 c open(46,file='whole.scr',form='unformatted') read(46)na1 do 1111 ia=1,na1 read(46)it,xc,yc,zc qc=ct(it) qi(ia)=qc xi(ia)=xc yi(ia)=yc zi(ia)=zc iti(ia)=it do 1111 ja=1,NAT r2=(xc-x(ja))**2+(yc-y(ja))**2+(zc-z(ja))**2 1111 cti(nc+ja)=cti(nc+ja)+qc/dsqrt(r2) close(46) write(3,*)' MD ave pot. ave positions/tinker charge pot' if(LIN)write(6,*) 1' MD ave pot. ave positions/tinker charge pot' do 1112 ja=1,NAT if(LIN)write(6,609)ja,izq(ja),c(nc+ja),cti(nc+ja) 1112 write(3,609)ja,izq(ja),c(nc+ja),cti(nc+ja) 609 format(2i4,2f15.6) c ***************************************************************** else c select NC charges for mapping in the first geometry ic=0 open(2,file='NAMES.LST',status='OLD') read(2,300)fn open(21,file=fn,status='OLD') call ratm(N1,N2,x,y,z,q,ct,izq,izt,nc0,LSX) rewind 21 read(21,*)natmd do 61 ia=1,natmd read(21,2121)s1,xa,ya,za,it if(it.eq.0.and.LSX)then if(s1.eq.'O')it=649 if(s1.eq.'H')it=650 c N,C,glycin types in amber.prm: if(s1.eq.'C')it=2 if(s1.eq.'N')it=1 izt(1)=7 izt(2)=6 izt(649)=8 izt(650)=1 endif if(ia.lt.N1.or.ia.gt.N2)then qc=ct(it) ic=ic+1 xi(ic)=xa yi(ic)=ya zi(ic)=za di(ic)=dmin(NAT,x,y,z,LPC,xa,ya,za,dx2,dy2,dz2,DPCX, 1DPCY,DPCZ) qi(ic)=qc endif 61 continue close(21) close(2) if(ic.lt.NC)call report('Not enough charges was found') c c order charges do 28 i=1,ic do 28 j=i+1,ic 28 if(di(j).lt.di(i))call sw(i,j,xi,yi,zi,qi,di) c if(LWR)write(3,*)'charge distance from molecule' if(LWR.and.LIN)write(6,*)'charge distance from molecule' do 26 i=1,NC if(LWR)write(3,2626)i,di(i)**0.5d0,qi(i) 26 if(LWR.and.LIN)write(6,2626)i,di(i)**0.5d0,qi(i) 2626 format(i6,2f12.6) endif endif c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c if(LAM.and.LAT)then write(3,*)' Tinker charges left for average mol. positions' if(LIN)write(6,*) 1 ' Tinker charges left for average mol. positions' do 120 i=1,NC 120 b(i)=qi(i) else if(LQ0)then write(3,*)'charges close to original seeked' do 4 i=1,NC 4 c(i)=qi(i) else write(3,*)'minimal charges seeked' endif c c set up the matrix do 14 i=1,M do 14 j=1,M 14 a(i,j)=0.0d0 do 15 i=1,NC a(i,i)=1.0d0 a(i,M)=1.0d0 15 a(M,i)=1.0d0 do 16 ia=1,NAT xa=x(ia) ya=y(ia) za=z(ia) do 16 ic=1,NC rij=dsqrt((xi(ic)-xa)**2+(yi(ic)-ya)**2+(zi(ic)-za)**2) a(NC+ia,ic)=1.0d0/rij 16 a(ic,NC+ia)=1.0d0/rij do 17 i=1,M do 17 j=1,M 17 a0(i,j)=a(i,j) if(LWR)call wm(a,M0,M) M1=M M2=M call dgetrf(M1,M2,a,M0,ipiv,info) if(info.eq.0)then if(LIN)write(6,*)' Pivoting OK' write(3,*)' 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(3,*)'successful inversion' if(info.eq.0.and.LIN)write(6,*)'successful inversion' if(info.lt.0)then write(3,*)info if(LIN)write(6,*)info call report('invalid argument') endif c tol=toi 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)*a0(ii,j) sum2=sum2+a(ii,i)*a0(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(3,*)'ic = ',ic write(3,*)i,j,sum,sum2 write(3,*)i,j,sum call report('Inversion error') endif enddo enddo c write(3,*)'Inversion Check OK' if(LIN)write(6,*)'Inversion Check OK' do 18 i=1,M sum=0.0d0 do 19 j=1,M 19 sum=sum+a(i,j)*c(j) 18 b(i)=sum endif write(3,88) if(LIN)write(6,88) 88 format(' charge old new') do 20 i=1,NC if(LIN)write(6,89)i,qi(i),b(i) 20 write(3,89)i,qi(i),b(i) 89 format(i6,2f10.2) open(4,file='GEO.X') write(4,*)' Atoms and charges from genchar' write(4,*)NAT+NC do 21 ia=1,NAT 21 write(4,3121)izq(ia),x(ia),y(ia),z(ia),q(ia) 3121 format(i6,3f12.6,' 0 0 0 0 0 0 0 ',f10.2) do 22 ic=1,NC iy=9 if(LAM)then iy=1 if(iti(ic).eq.649)iy=8 endif 22 write(4,3121)iy,xi(ic),yi(ic),zi(ic),b(ic) close(4) write(3,*)' GEO.X written' if(LIN)write(6,*)' GEO.X written' open(4,file='GEO.INP') write(4,400) 400 format('%mem=24000000',/, 1'%nproc=1',/, 2'%rwf=1,249MW,2,249MW,3,249MW,4,249MW,5,249MW,6,-1',/, 3'#BPW91/6-31G** force nosymm fmm charge',/, 4/, 5' Atoms and charges from genchar',/, 6/, 7'0 1 ') do 23 ia=1,NAT 23 write(4,2122)izq(ia),x(ia),y(ia),z(ia) 2122 format(i6,3f12.6) write(4,*) do 24 ic=1,NC 24 write(4,2123)xi(ic),yi(ic),zi(ic),b(ic) 2123 format(4f12.6) write(4,*) close(4) write(3,*)' GEO.INP written' if(LIN)write(6,*)' GEO.INP written' write(3,*)'test pot old new' if(LIN)write(6,*)'test pot old new' do 25 ia=1,nat pot=0.0d0 do 27 ic=1,nc 27 pot=pot+b(ic)/dsqrt((x(ia)-xi(ic))**2+(y(ia)-yi(ic))**2 1 +(z(ia)-zi(ic))**2) if(LIN)write(6,2525)ia,c(ia+NC),pot 25 write(3,2525)ia,c(ia+NC),pot 2525 format(i6,2f12.6) close(3) stop end c =========================================== subroutine iopar(LAV,LAM,LAT,LPC,LSP,LQ0,LWR,LIN,NMC,NC, 1DPCX,DPCY,DPCZ,NAK,NAM,itk,ntk,ntm,N1,N2,fnpar, 2maxit,tol,toi,LSX) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) character*3 s3 character*80 fnpar logical LPC,LWR,LAV,LAM,LAT,LIN,LQ0,LSP,LSX dimension itk(*),ntk(*) toi=1.0d-5 open(2,file='CHAR.OPT',status='OLD') 1 read(2,200,ERR=2,END=2)s3 200 format(a3) if(s3.eq.'LAV')read(2,*)LAV if(s3.eq.'LAM')read(2,*)LAM if(s3.eq.'LAT')read(2,*)LAT if(s3.eq.'LPC')read(2,*)LPC if(s3.eq.'LSP')read(2,*)LSP if(s3.eq.'LSX')read(2,*)LSX if(s3.eq.'LQ0')read(2,*)LQ0 if(s3.eq.'LWR')read(2,*)LWR if(s3.eq.'INT')read(2,*)LIN if(s3.eq.'NMC')read(2,*)NMC if(s3.eq.'NCH')read(2,*)NC if(s3.eq.'DPC')read(2,*)DPCX,DPCY,DPCZ if(s3.eq.'NAK')read(2,*)NAK if(s3.eq.'NAM')read(2,*)NAM if(s3.eq.'ITK')read(2,*)iitk,itk(iitk) if(s3.eq.'NTK')read(2,*)intk,ntk(intk) if(s3.eq.'NTM')read(2,*)ntm if(s3.eq.'N1N')read(2,*)N1,N2 if(s3.eq.'PAR')read(2,'(a)')fnpar if(s3.eq.'MAX')read(2,*)maxit if(s3.eq.'TOL')read(2,*)tol if(s3.eq.'TOI')read(2,*)toi if(s3.eq.'END')goto 2 goto 1 2 close(2) write(3,*)'Options read from CHAR.OPT' if(LIN)write(6,*)'Options read from CHAR.OPT' return end c =========================================== subroutine makegr(NCO,NO,tkx,tky,tkz,iatom,maxit,tol) parameter (natall=500000,nt0=4) real*4 rijx(natall),rijy(natall),rijz(natall),tkx(*), 1tky(*),tkz(*),s0,tol4 character*8 nte character*80 fnpar real*8 tol,DPCX,DPCY,DPCZ,toi integer*4 NCO,NO,ia,i,j,ii,iatom(*),maxit,it common/bigc/rijx,rijy,rijz logical LPC,LWR,LAV,LAM,LAT,LIN,LQ0,LSP,LSX dimension itk(nt0),ntk(nt0) common/inter/LIN c c initialize groups tol4=real(tol) RNCO=real(NCO) do 1 ia=1,NO*NCO 1 iatom(ia)=ia s0= si(NO,iatom,NCO,rijx,rijy,rijz,tkx,tky,tkz) write(3,*)'Initial dispersion = ',s0 if(LIN)write(6,*)'Initial dispersion = ',s0 it0=time() it=0 c 44 sold=si(NO,iatom,NCO,rijx,rijy,rijz,tkx,tky,tkz) it=it+1 call iopar(LAV,LAM,LAT,LPC,LSP,LQ0,LWR,LIN,NMC,NC, 1DPCX,DPCY,DPCZ,NAK,NAM,itk,ntk,ntm,N1,N2,fnpar, 2maxit,tol,toi,LSX) tol4=real(tol) do 4 i1=1,NO t1x=tkx(i1) t1y=tky(i1) t1z=tkz(i1) ii1=NCO*(i1-1) do 4 j1=1,NCO ia1=iatom(ii1+j1) rijx1=rijx(ia1) rijy1=rijy(ia1) rijz1=rijz(ia1) do 4 i2=i1+1,NO t2x=tkx(i2)-t1x t2y=tky(i2)-t1y t2z=tkz(i2)-t1z ii2=NCO*(i2-1) do 4 j2=1,NCO ia2=iatom(ii2+j2) rijx2=rijx(ia2)-rijx1 rijy2=rijy(ia2)-rijy1 rijz2=rijz(ia2)-rijz1 rt=rijx2*t2x+rijy2*t2y+rijz2*t2z rh=(rijx2**2+rijy2**2+rijz2**2)/RNCO c c look if switching atoms 1 2 helps if(rt.lt.rh)then iatom(ii2+j2)=ia1 iatom(ii1+j1)=ia2 ia1=ia2 ia2=iatom(ii2+j2) rijx1=rijx(ia1) rijy1=rijy(ia1) rijz1=rijz(ia1) tkx(i1)=tkx(i1)+rijx2/RNCO tky(i1)=tky(i1)+rijy2/RNCO tkz(i1)=tkz(i1)+rijz2/RNCO tkx(i2)=tkx(i2)-rijx2/RNCO tky(i2)=tky(i2)-rijy2/RNCO tkz(i2)=tkz(i2)-rijz2/RNCO t1x=tkx(i1) t1y=tky(i1) t1z=tkz(i1) t2x=tkx(i2)-t1x t2y=tky(i2)-t1y t2z=tkz(i2)-t1z endif 4 continue s=si(NO,iatom,NCO,rijx,rijy,rijz,tkx,tky,tkz) write(3,6500)it,s,sold-s,tol4 if(LIN)write(6,6500)it,s,sold-s,tol4 6500 format(i6,' s = ',f15.6,', ds = ',f15.6,', tol = ',f10.4) if(sold-s.le.tol4)then write(3,*)' Change below limit achieved' if(LIN)write(6,*)' Change below limit achieved' goto 45 endif if(it.ge.maxit)then write(3,*)' Maximum number of iterations overflow' if(LIN)write(6,*)' Maximum number of iterations overflow' goto 45 endif goto 44 45 write(3,*)' time ',time()-it0 if(LIN)write(6,*)' time ',time()-it0 open(8,file='scr') do 8 i=1,NO rewind 8 write(8,8000)i,'.x' 8000 format(i6,a2) rewind 8 read(8,8001)nte 8001 format(a8) do 101 ist=1,len(nte) 101 if(nte(ist:ist).ne.' ')goto 102 102 open(45,file=nte(ist:len(nte))) write(45,8001)nte write(45,*)NCO+1 ii=NCO*(i-1) do 9 j=1,NCO ia=iatom(ii+j) 9 write(45,4545)rijx(ia),rijy(ia),rijz(ia) 4545 format(' 1 ',3f12.6) write(45,4546)tkx(i),tky(i),tkz(i) 4546 format(' 8 ',3f12.6) 8 close(45) close(8) return end c =========================================== function si(NO,iatom,NCO,rijx,rijy,rijz,tkx,tky,tkz) integer*4 NO,NCO,iatom(*) dimension rijx(*),rijy(*),rijz(*),tkx(*),tky(*), 1tkz(*) RNCO=real(NCO) do 111 i=1,NO sx=0.0 sy=0.0 sz=0.0 ii=NCO*(i-1) do 211 j=1,NCO ia=iatom(ii+j) sx=sx+rijx(ia) sy=sy+rijy(ia) 211 sz=sz+rijz(ia) tkx(i)=sx/RNCO tky(i)=sy/RNCO 111 tkz(i)=sz/RNCO c s=0.0 do 51 i=1,NO tx=tkx(i) ty=tky(i) tz=tkz(i) ii=NCO*(i-1) do 51 j=1,NCO ia=iatom(ii+j) 51 s=s+(tx-rijx(ia))**2+(ty-rijy(ia))**2+(tz-rijz(ia))**2 si=s return end c =========================================== SUBROUTINE mopti IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) COMMON/MATVAR/A(3,3),RTOL,ITER parameter (NAM0=10,natall=500000) dimension xo(NAM0),yo(NAM0),zo(NAM0),xn(NAM0),yn(NAM0),zn(NAM0), 1xh(NAM0),yh(NAM0),zh(NAM0),iatom(natall), 2tlistx(natall,NAM0),tlisty(natall,NAM0),tlistz(natall,NAM0), 3itlist(natall,NAM0) common/mv2/xo,yo,zo,xh,yh,zh,tlistx,tlisty,tlistz,xn,yn,zn, 1NAM,NMC,ii1,iatom,itlist logical LIN common/inter/LIN CALL DOMATRIX(IERR) if(IERR.ne.0)call report('matrix not found') if(LIN)write(6,*)ITER write(3,*)ITER do 1 ia=1,NAM xn(ia)=A(1,1)*xo(ia)+A(1,2)*yo(ia)+A(1,3)*zo(ia) yn(ia)=A(2,1)*xo(ia)+A(2,2)*yo(ia)+A(2,3)*zo(ia) 1 zn(ia)=A(3,1)*xo(ia)+A(3,2)*yo(ia)+A(3,3)*zo(ia) 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) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3) COMMON/MATVAR/A(3,3),RTOL,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=0.0000001d0 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.0.00000000001d0)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) COMMON/MATVAR/A(3,3),RTOL,ITER parameter (NAM0=10,natall=500000) dimension xo(NAM0),yo(NAM0),zo(NAM0),xn(NAM0),yn(NAM0),zn(NAM0), 1xh(NAM0),yh(NAM0),zh(NAM0),iatom(natall), 2tlistx(natall,NAM0),tlisty(natall,NAM0),tlistz(natall,NAM0), 3itlist(natall,NAM0) common/mv2/xo,yo,zo,xh,yh,zh,tlistx,tlisty,tlistz,xn,yn,zn, 1NAM,NMC,ii1,iatom,itlist 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 c c rotate the molecule DO 1 I=1,NAM xn(I)=A(1,1)*xo(I)+A(1,2)*yo(I)+A(1,3)*zo(I) yn(I)=A(2,1)*xo(I)+A(2,2)*yo(I)+A(2,3)*zo(I) 1 zn(I)=A(3,1)*xo(I)+A(3,2)*yo(I)+A(3,3)*zo(I) c RNAM=dble(NAM) c loop over all configurations in the group do 2 j=1,NMC iig=iatom(ii1+j-1) tx=0.0d0 ty=0.0d0 tz=0.0d0 do 3 ia=1,NAM xh(ia)=tlistx(iig,ia) yh(ia)=tlisty(iig,ia) zh(ia)=tlistz(iig,ia) tx=tx+xh(ia) ty=ty+yh(ia) 3 tz=tz+zh(ia) tx=tx/RNAM ty=ty/RNAM tz=tz/RNAM do 4 ia=1,NAM xh(ia)=xh(ia)-tx yh(ia)=yh(ia)-ty 4 zh(ia)=zh(ia)-tz do 2 ia=1,NAM 2 R=R+(xh(ia)-xn(ia))**2+(yh(ia)-yn(ia))**2+(zh(ia)-zn(ia))**2 FU=R RETURN END C c c what is the minimum square distance from the molecule: function dmin(NAT,x,y,z,LPC,xa,ya,za,dx2,dy2,dz2,DPCX, 1DPCY,DPCZ) real*8 dmin,x(*),y(*),z(*),xj,yj,zj,rc,dmm, 1xa,ya,za,dx2,dy2,dz2,DPCY,DPCZ,DPCX integer*4 ja,NAT logical LPC dmm=1.0d20 do 101 ja=1,NAT xj=x(ja) yj=y(ja) zj=z(ja) if(LPC)then if(xa-xj.gt.dx2)xa=xa-DPCX if(xj-xa.gt.dx2)xa=xa+DPCX if(ya-yj.gt.dy2)ya=ya-DPCY if(yj-ya.gt.dy2)ya=ya+DPCY if(za-zj.gt.dz2)za=za-DPCZ if(zj-za.gt.dz2)za=za+DPCZ endif rc=(xj-xa)**2+(yj-ya)**2+(zj-za)**2 101 if(rc.lt.dmm)dmm=rc dmin=dmm return end c =========================================== subroutine ratm(N1,N2,x,y,z,q,ct,izq,izt,nc0,LSX) integer*4 N1,N2,ia,it,l,izq(*),izt(*),nc0 real*8 xa,ya,za,x(*),y(*),z(*),q(*),ct(*) character*1 s1 logical LSX c read(21,*) do 51 ia=1,N2 read(21,2121)s1,xa,ya,za,it 2121 format(6x,2x,a1,2x,3f12.6,i6) if(it.eq.0.and.LSX)then if(s1.eq.'O')it=649 if(s1.eq.'H')it=650 c N,C,glycin types in amber.prm: if(s1.eq.'C')it=2 if(s1.eq.'N')it=1 izt(1)=7 izt(2)=6 izt(649)=8 izt(650)=1 endif if(ia.ge.N1)then l=ia-N1+1 x(l)=xa y(l)=ya z(l)=za q(l)=ct(it) if(it.lt.1.or.it.gt.nc0)then write(3,*)ia,it,nc0 call report('Invalid it.') endif izq(l)=izt(it) endif 51 continue 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 subroutine report(s) character*(*) s logical LIN common/inter/LIN write(3,*)s if(LIN)write(6,*)s stop end subroutine co(NMC,NC,LPC,DPCX,DPCY,DPCZ,N1,N2,fnpar,NAT,M) implicit real*8 (a-h,o-z) implicit integer*4 (i-n) logical LPC character*80 fnpar character*3 yy logical LIN common/inter/LIN yy=' NO' if(LPC)yy='YES' write(3,6000)NMC,NC,yy,DPCX,DPCY,DPCZ,N1,N2,fnpar(1:40),NAT,M if(LIN) 1write(6,6000)NMC,NC,yy,DPCX,DPCY,DPCZ,N1,N2,fnpar(1:40),NAT,M 6000 format(' Number of MD configurations: ',I6,/, 1 ' Number of charges : ',I6,/, 2 ' Periodic bondary conditions: ',3x,A3,/, 3 ' Periodic bondary X : ',f12.6,/, 4 ' Periodic bondary y : ',f12.6,/, 5 ' Periodic bondary z : ',f12.6,/, 6 ' Solvated molecule within : ',i6,' - ',i6,/, 7 ' MD force field : ',a40,/, 8 ' Number of atoms/pot. sites : ',i6,/, 9 ' Dimension of inv. matrix : ',i6,/) return end subroutine wm(a,M0,M) integer*4 M0,M,i,j real*8 a(M0,M0) logical LIN common/inter/LIN do 1 i=1,M write(3,600)i if(LIN)write(6,600)i 600 format(' i = ',i6) if(LIN)write(6,606)(a(i,j),j=1,M) 1 write(3,606)(a(i,j),j=1,M) 606 format(6f12.6) return end subroutine swm(i,j,xi,yi,zi,qi,di,alistx,alisty,alistz,NAM0, 1NAM,listi,nc00) integer*4 i,j,NAM0,NAM,k,listi(nc00,NAM0) real*8 xi(*),yi(*),zi(*),qi(*),di(*), 1alistx(nc00,NAM0),alisty(nc00,NAM0),alistz(nc00,NAM0) call swr(xi(i),xi(j)) call swr(yi(i),yi(j)) call swr(zi(i),zi(j)) call swr(di(i),di(j)) call swr(qi(i),qi(j)) do 1 k=1,NAM call swi( listi(i,k), listi(j,k)) call swr(alistx(i,k),alistx(j,k)) call swr(alisty(i,k),alisty(j,k)) 1 call swr(alistz(i,k),alistz(j,k)) return end subroutine sw(i,j,xi,yi,zi,qi,di) integer*4 i,j real*8 xi(*),yi(*),zi(*),qi(*),di(*) call swr(xi(i),xi(j)) call swr(yi(i),yi(j)) call swr(zi(i),zi(j)) call swr(di(i),di(j)) call swr(qi(i),qi(j)) return end subroutine swr(a,b) real*8 a,b,c c=a a=b b=c return end subroutine swi(a,b) integer*4 a,b,c c=a a=b b=c return end