PROGRAM FITNMA C finds fitting coefficients for solvation effect c on force field belonging to atoms 1...N4 c based on potential on atoms 1..N6 c nat0 - maximum number of atoms c m0 - max number of clusters c nd0 - max dimension of the inversion problem IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (nat0=100,MX3=3*nat0,MX=nat0,m0=100,nd0=10) CHARACTER*80 sname,snamel(m0),xnamel(m0),tnamel(m0),pnamel(m0) DIMENSION ZIM(MX),ZM(MX3),FCAR(MX3,MX3), 2SMAT(MX3,MX3),Q(MX3),fcars(MX3,MX3),fi(m0,nat0),r(nat0,3), 3ity(nat0),nal(m0),i6(m0,nat0),i4(m0,nat0),qq(m0,MX3), 4a(nd0,nd0),b(nat0,nd0),c(nd0),am(nd0,nd0),ipiv(nd0), 5work(3*nd0*nd0),E1(MX3),P1(nat0,3,3),fi_ref(nat0), 6P(nat0,3,3),AAT(nat0,3,3),qqx(m0,nat0,3,3),qpx(m0,nat0,3,3,3), 7bx(nat0,3,nd0,3),e0(3,3),bp(nat0,3,nd0,3,3), 8alpha(MX3,3,3),AT(MX3,3,3,3),G(MX3,3,3),qt(MX3,3,3) LOGICAL lfdp(3) DATA IERR/0/ DATA qh,qo/0.42d0,-0.84d0/ C WRITE(6,60000) 60000 FORMAT(' FORCE FIELD SOLVENT PERTURBATION',/,/, 1 ' Input : FITNMA.INP ',/, 1 ' Output: BIJ.SCR ',/) C OPEN(16,FILE='FITNMA.INP',status='old') read(16,1600)sname write(6,1600)sname 1600 format(a80) read(16,*)m write(6,1601)m 1601 format(i5,' md configurations') if(m.gt.m0)call report(' Too many configurations') read(16,*)n4,n6 write(6,1602)n4,n6 1602 format(i5,' perturbed atoms, ',i2,' potential sites') if(n6.gt.nd0-1)call report('Too many pot. sites') c c read masses, just for units of the FCs: read(16,*)(zim(i),i=1,n4) c c read whether to fir fc, P, A: read(16,*)(lfdp(i),i=1,3) if(lfdp(1))then write(6,*)'FC requested' else write(6,*)'FC ommitted' endif if(lfdp(2))then write(6,*)'P requested' else write(6,*)'P ommitted' endif if(lfdp(3))then write(6,*)'A requested' else write(6,*)'A ommitted' endif DO 220 i=1,n4 DO 220 j=1,3 220 ZM(3*(i-1)+j)=1.0d0/DSQRT(zim(i)) do 1 i=1,m c FC, dipole derivaties, polarizability derivatives, coords: if(lfdp(1))read(16,1600)snamel(i) if(lfdp(2))read(16,1600)tnamel(i) if(lfdp(3))read(16,1600)pnamel(i) read(16,1600)xnamel(i) read(16,*)nal(i),(i6(i,j),j=1,n6),(i4(i,j),j=1,n4) 1 if(nal(i).gt.nat0)call report('Too many atoms') close(16) c ========================================================= do 2 i=1,m write(6,1600)xnamel(i) c c read geometry OPEN(20,FILE=xnamel(i),STATUS='OLD') read(20,*) read(20,*)nat do 4 ia=1,nat 4 read(20,*)ity(ia),(r(ia,ix),ix=1,3) close(20) c c local coordinate system do 50 ix=1,3 e0(1,ix)=r(i4(i,4),ix)-r(i4(i,3),ix) 50 e0(3,ix)=r(i4(i,2),ix)-r(i4(i,3),ix) call norm(e0,1) e0(2,1)=e0(1,2)*e0(3,3)-e0(1,3)*e0(3,2) e0(2,2)=e0(1,3)*e0(3,1)-e0(1,1)*e0(3,3) e0(2,3)=e0(1,1)*e0(3,2)-e0(1,2)*e0(3,1) call norm(e0,2) e0(3,1)=e0(1,2)*e0(2,3)-e0(1,3)*e0(2,2) e0(3,2)=e0(1,3)*e0(2,1)-e0(1,1)*e0(2,3) e0(3,3)=e0(1,1)*e0(2,2)-e0(1,2)*e0(2,1) c c calculate potentials on atoms 1..N6 do 5 ii=1,N6 x=r(i6(i,ii),1) y=r(i6(i,ii),2) z=r(i6(i,ii),3) sum=0.0d0 do 6 jj=1,nat if(ity(jj).eq.8)then xo=r(jj,1) yo=r(jj,2) zo=r(jj,3) nh=0 do 7 kk=1,nat if(ity(kk).eq.1)then xh=r(kk,1) yh=r(kk,2) zh=r(kk,3) d2=(xo-xh)**2+(yo-yh)**2+(zo-zh)**2 if(d2.gt.0.64d0.and.d2.lt.1.44d0)then nh=nh+1 if(nh.gt.2)call report('Too many Hs around one O') if(nh.eq.1)then x1=xh y1=yh z1=zh else x2=xh y2=yh z2=zh endif endif endif 7 continue if(nh.eq.2)then ro=dsqrt((xo-x)**2+(yo-y)**2+(zo-z)**2) r1=dsqrt((x1-x)**2+(y1-y)**2+(z1-z)**2) r2=dsqrt((x2-x)**2+(y2-y)**2+(z2-z)**2) sum=sum+qh/r1+qh/r2+qo/ro endif endif 6 continue 5 fi(i,ii)=sum write(6,*)' Potentials:' write(6,1607)(fi(i,ii),ii=1,N6) 1607 format(8f10.3) if(i.eq.1)then do ii=1,N6 fi_ref(ii)=fi(i,ii) enddo endif write(6,*)' Potentials - difference vs reference (first cluster):' do ii=1,N6 fi(i,ii)=fi(i,ii)-fi_ref(ii) enddo write(6,1607)(fi(i,ii),ii=1,N6) c C Read-in cartesian FF if(lfdp(1))then OPEN(20,FILE=snamel(i),STATUS='OLD') CALL READFF(MX3,3*nal(i),FCAR) CLOSE(20) c rewrite into small matrix do 3 ii=1,n4 ia=i4(i,ii) do 3 ix=1,3 iin=3*(ii-1)+ix i1 =3*(ia-1)+ix do 3 jj=1,n4 ja=i4(i,jj) do 3 jx=1,3 jin=3*(jj-1)+jx j1 =3*(ja-1)+jx 3 fcars(iin,jin)=FCAR(i1,j1) c calculate eigenvalues IMOD=2 NA=3*N4 DO 214 i1=1,NA DO 214 j1=1,NA 214 fcars(i1,j1)=fcars(i1,j1)*zm(j1)*zm(i1) CALL TRED12(MX3,NA,fcars,SMAT,Q,IMOD,IERR,E1) IF(IERR.NE.0)call report(' Diagonalization error') DO 310 Iz=1,NA t=Q(Iz) qq(i,Iz)=t Q(Iz)=SQRT(ABS(t))*1302.828d0 310 if(t.lt.0.0d0)Q(Iz)=-Q(Iz) write(6,*)' Eigenvalues:' write(6,1603)(Q(II),II=1,NA) 1603 format(8f10.1) endif c c read-in dipole derivatives if(lfdp(2))then OPEN(15,FILE=tnamel(i),STATUS='OLD') call READTEN(nat0,P,AAT,NAT) close(15) c c transform dipole derivatives to local system write(6,*)'dipole der in local system:' do 51 ii=1,N4 ia=i4(i,ii) c c first index do 511 k=1,3 do 511 j=1,3 P1(ia,j,k)=0.0d0 do 511 l=1,3 511 P1(ia,j,k)=P1(ia,j,k)+e0(j,l)*P(ia,l,k) c c second index do 512 k=1,3 do 512 j=1,3 P(ia,j,k)=0.0d0 do 512 l=1,3 512 P(ia,j,k)=P(ia,j,k)+e0(k,l)*P1(ia,j,l) do 51 k=1,3 do 51 l=1,3 write(6,643)ii,k,l,P(ia,k,l) 643 format(3i3,f12.3) 51 qqx(i,ii,k,l)=P(ia,k,l) endif c c read-in polarizability derivatives if(lfdp(3))then OPEN(2,FILE=pnamel(i),STATUS='OLD') call READTTT(MX3,NAT,alpha,G,AT) close(2) c c transform to local system write(6,*)'pol der in local system:' do 52 ii=1,N4 ia=i4(i,ii) c c first(atomic) index do 611 i1=1,3 i0=i1+3*(ia-1) do 611 i2=1,3 do 611 i3=1,3 qt(i0,i2,i3)=0.0d0 do 611 l=1,3 611 qt(i0,i2,i3)=qt(i0,i2,i3)+e0(i1,l)*alpha(l+3*(ia-1),i2,i3) c c second index do 612 i1=1,3 i0=i1+3*(ia-1) do 612 i2=1,3 do 612 i3=1,3 alpha(i0,i2,i3)=0.0d0 do 612 l=1,3 612 alpha(i0,i2,i3)=alpha(i0,i2,i3)+e0(i2,l)*qt(i0,l,i3) c c third index do 613 i1=1,3 i0=i1+3*(ia-1) do 613 i2=1,3 do 613 i3=1,3 qt(i0,i2,i3)=0.0d0 do 613 l=1,3 613 qt(i0,i2,i3)=qt(i0,i2,i3)+e0(i3,l)*alpha(i0,i2,l) c do 52 mm=1,3 do 52 k=1,3 do 52 l=1,3 write(6,644)ii,mm,k,l,qt(mm+3*(ia-1),k,l) 644 format(4i3,f12.3) 52 qpx(i,ii,mm,k,l)=qt(mm+3*(ia-1),k,l) endif c 2 continue c ========================================================= c c ========================================================= c c set-up matrix (A B = C, B=A^-1 C) nd=n6+1 do 10 i1=1,n6 do 10 j1=1,n6 sum=0.0d0 do 11 i=1,m 11 sum=sum+fi(i,i1)*fi(i,j1) 10 a(i1,j1)=sum do 12 i1=1,n6 a(i1,nd)=1.0d0 12 a(nd,i1)=1.0d0 a(nd,nd)=0.0d0 c c invert call inv(a,nd,nd0,ipiv,work,am) c if(lfdp(1))then c c form c for each eigenvalue,FF FIT do 13 i=1,3*n4 write(6,*)'Eigenvalue ',i do 14 j=1,n6 sum=0.0d0 do 15 jj=1,m 15 sum=sum+(qq(jj,i)-qq(1,i))*fi(jj,j) 14 c(j)=sum c(nd)=0.0d0 c c get B do 16 ii=1,n6 sum=0.0d0 do 17 jj=1,nd 17 sum=sum+am(ii,jj)*c(jj) 16 b(i,ii)=sum 13 write(6,1604)(b(i,ii),ii=1,n6) 1604 format(6f12.6) c c check fit write(6,*)' Old and fitted constants:' do 18 i=1,m do 18 j=1,7 old=qq(i,j)-qq(1,j) sum=0.0d0 do 19 ii=1,n6 19 sum=sum+fi(i,ii)*b(j,ii) 18 write(6,1605)i,j,qq(i,j),old,sum 1605 format(2i3,3f12.6) endif c if(lfdp(2))then c c form c for each eigenvalue and dipole component, dipole FIT do 131 i=1,n4 do 131 iy=1,3 do 131 ix=1,3 write(6,*)'At ',i,' coord ',iy,' dip ',ix do 141 j=1,n6 sum=0.0d0 do 151 jj=1,m 151 sum=sum+(qqx(jj,i,iy,ix)-qqx(1,i,iy,ix))*fi(jj,j) 141 c(j)=sum c(nd)=0.0d0 c c get Bx do 161 ii=1,n6 sum=0.0d0 do 171 jj=1,nd 171 sum=sum+am(ii,jj)*c(jj) 161 bx(i,iy,ii,ix)=sum 131 write(6,1604)(bx(i,iy,ii,ix),ii=1,n6) write(6,*)' Old and fitted dipole constants:' do 181 i=1,m do 181 j=1,4 do 181 jx=1,3 do 181 ix=1,3 old=qqx(i,j,jx,ix)-qqx(1,j,jx,ix) sum=0.0d0 do 191 ii=1,n6 191 sum=sum+fi(i,ii)*bx(j,jx,ii,ix) 181 write(6,16051)i,j,jx,ix,qqx(i,j,jx,ix),old,sum 16051 format(4i3,3f12.5) endif c if(lfdp(3))then c c form c for each component, polarizability FIT do 132 i=1,n4 do 132 iy=1,3 do 132 ix=1,3 do 132 iz=1,3 write(6,*)'At ',i,' coord ',iy,' pol ',ix,iz do 142 j=1,n6 sum=0.0d0 do 152 jj=1,m 152 sum=sum+(qpx(jj,i,iy,ix,iz)-qpx(1,i,iy,ix,iz))*fi(jj,j) 142 c(j)=sum c(nd)=0.0d0 c c get Bp do 162 ii=1,n6 sum=0.0d0 do 172 jj=1,nd 172 sum=sum+am(ii,jj)*c(jj) 162 bp(i,iy,ii,ix,iz)=sum 132 write(6,1604)(bp(i,iy,ii,ix,iz),ii=1,n6) write(6,*)' Old and fitted polarizability' do 182 i=1,m do 182 j=1,4 do 182 jx=1,3 do 182 ix=1,3 do 182 iz=1,3 old=qpx(i,j,jx,ix,iz)-qpx(1,j,jx,ix,iz) sum=0.0d0 do 192 ii=1,n6 192 sum=sum+fi(i,ii)*bp(j,jx,ii,ix,iz) 182 write(6,16052)i,j,jx,ix,iz,qpx(i,j,jx,ix,iz),old,sum 16052 format(5i3,3f12.5) endif c c save to disk open(55,file='BIJ.SCR',form='unformatted') write(55)n4,n6,(zim(i),i=1,n4),(lfdp(i),i=1,3) if(lfdp(1))write(55)((b(i,j),i=1,3*n4),j=1,n6) if(lfdp(2))write(55) 1((((bx(i,i3,j,ix),i=1,n4),i3=1,3),j=1,n6),ix=1,3) if(lfdp(3))write(55) 1(((((bp(i,i3,j,ix,iz),i=1,n4),i3=1,3),j=1,n6),ix=1,3),iz=1,3) write(55)n6,(fi_ref(i),i=1,n6) close(55) open(55,file='BIJ.TXT') write(55,555)n4,n6 write(55,552)(zim(i),i=1,n4) write(55,553)(lfdp(i),i=1,3) 555 format(2i6) 552 format(8f10.1) 553 format(3l2) if(lfdp(1))then do i=1,3*n4 write(55,554)(b(i,j),j=1,n6) 554 format(8G10.2) enddo endif close(55) c call report('Fitting constants have been saved to BIJ.SCR.') end subroutine norm(v,i) real*8 v(3,3),an integer*4 i an=dsqrt(v(i,1)**2+v(i,2)**2+v(i,3)**2) v(i,1)=v(i,1)/an v(i,2)=v(i,2)/an v(i,3)=v(i,3)/an return end subroutine report(s) character*(*) s write(6,*)s stop end c SUBROUTINE READTTT(N30,NAT,ALPHA,G,A) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ALPHA(N30,3,3),G(N30,3,3),A(N30,3,3,3) READ(2,*) READ(2,*)NAT READ(2,*) READ(2,*) DO 1 I=1,3 READ(2,*) DO 1 L=1,NAT DO 1 IX=1,3 IND=3*(L-1)+IX 1 READ(2,*)LDUM,IXDUM,(ALPHA(IND,I,J),J=1,3) READ(2,*) READ(2,*) DO 2 I=1,3 READ(2,*) DO 2 L=1,NAT DO 2 IX=1,3 IND=3*(L-1)+IX c last index magnetic, inner 2 READ(2,*)LDUM,IXDUM,(G(IND,I,J),J=1,3) READ(2,*) READ(2,*) DO 3 I=1,3 DO 3 J=1,3 READ(2,*) DO 3 L=1,NAT DO 3 IX=1,3 IND=3*(L-1)+IX 3 READ(2,*)LDUM,IXDUM,(A(IND,I,J,K),K=1,3) RETURN END c SUBROUTINE READTEN(N0,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N0,3,3),A(N0,3,3) READ(15,*) DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) DO 220 L=1,NAT DO 220 J=1,3 220 READ (15,*) (A(L,J,I),I=1,3) DO 230 L=1,NAT DO 230 J=1,3 230 READ (15,*) DO 100 L=1,NAT DO 100 J=1,3 100 READ (15,*) RETURN END C SUBROUTINE TRED12(MX3,N,A,Z,D,IEIG,IERR,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(MX3,N),Z(MX3,N),D(MX3),E(MX3) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (MX3,N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(MX3,N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(MX3,N),D(MX3),E(MX3) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END C SUBROUTINE READFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX3,N) N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) C CONST=4.359828d0/0.5291772d0/0.5291772d0 DO 3 I=1,N DO 3 J=I,N 3 FCAR(J,I)=FCAR(J,I)*CONST DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) WRITE(*,*)' Cartesian FF read in ... ' RETURN END C c ====================================================================== subroutine inv(fx,M,M0,ipiv,work,a) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION fx(M0,M0),a(M0,M0),ipiv(M0),work(3*M0*M0) do 3 i=1,M do 3 j=i,M a(j,i)=fx(i,j) 3 a(i,j)=fx(i,j) M1=M M2=M call dgetrf(M1,M2,a,M0,ipiv,info) if(info.eq.0)then write(6,*)' Pivoting OK' else call report('pivoting failed') endif c lwork=3*M*M call dgetri(M,a,M0,ipiv,work,lwork,info) c if(info.eq.0)write(6,*)'successful inversion' if(info.lt.0)then call report('invalid argument') endif c tol=1.0d-7 nmo=M do i=1,nmo do j=1,nmo sum=0.0d0 sum2=0.0d0 do ii=1,nmo sum=sum+a(i,ii)*fx(ii,j) sum2=sum2+a(ii,i)*fx(j,ii) enddo ic=0 if(abs(sum ).gt.tol.and.i.ne.j)ic=1 if(abs(sum2).gt.tol.and.i.ne.j)ic=2 if(abs(sum -1.0d0).gt.tol.and.i.eq.j)ic=3 if(abs(sum2-1.0d0).gt.tol.and.i.eq.j)ic=4 if(ic.ne.0)then write(6,6990)((a(i1,j1),j1=1,M),i1=1,M) write(6,*)'ic = ',ic write(6,*)i,j,sum,sum2 write(6,6990)((fx(i1,j1),j1=1,M),i1=1,M) 6990 format(7f10.4) call report('Inversion error') endif enddo enddo c write(6,*)'Inversion Check OK' return end c 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