PROGRAM BUMATC6 c get the benzene stretching constants IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=435) COMMON/BMAT/IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,X(3,MX/3), 1B(MX,MX),NOB0,IS,UM(MX,MX) DIMENSION B0(MX,MX),LV(MX) logical lred call readopt(tol,lred,rf) WRITE(*,*)' (SYMMETRIZED) B-MATRIX GENERATION' WRITE(*,*)' FILE.UMA ... INPUT FILE' DO 1001 I=1,MX DO 1001 J=1,MX 1001 B(I,J)=0.0D0 OPEN(15,FILE='FILE.UMA') READ(15,*) READ(15,*) READ(15,*)NOAT,NOB0 IF(NOAT.GT.MX/3)THEN WRITE(*,*)' Too many atoms, redimension the program' CLOSE(15) STOP ENDIF WRITE(*,*)NOAT,NOB0 61000 format(/,' Warning: Number of Coordinates <> 3*NAT-6!',/) NA=NOAT*3 if(NOB0.ne.NA-6)write(6,61000) DO 20 I=1,NOAT READ(15,1020)(X(J,I),J=1,3) 1020 FORMAT(3F12.6) 20 CONTINUE IER=0 ISCAN=0 IS=0 DO 277 NOB=1,NOB0 READ(15,*)ICODE,N1,N2,N3,N4,N5,N6,FACTOR write(6,1050)ICODE,N1,N2,N3,N4,N5,N6,FACTOR 1050 FORMAT(7I4,F12.6) IF(DABS(FACTOR).LT.1.0D-6)FACTOR=1.0D0 IC=IABS(ICODE) IF(IC.EQ.0.OR.IC.GT.6) GO TO 220 GO TO (1,2,3,4,5,6),IC 220 WRITE(*,1070)NOB,NOB0,ICODE GO TO 200 1 CALL BOST GO TO 60 2 CALL BEND GO TO 60 3 CALL OPLA GO TO 60 4 CALL TORS GO TO 60 5 IS=IS+1 6 CALL LIBE 60 IF(IER.NE.0)GO TO 190 IF(DABS(FACTOR-1.0D0).LT.1.0D-8)GO TO 105 IF(IC.EQ.5)GO TO 80 70 ISW=0 I=NOB GO TO 90 80 ISW=1 I=NOB0+IS 90 DO 100 J=1,NA 100 B(I,J)=B(I,J)*FACTOR IF(ISW.EQ.1)GO TO 70 105 IF(ICODE.GT.0)GO TO 277 C IF(ICODE.EQ.-5)GO TO 110 C ISW=0 C I=NOB-1 C GO TO 130 C 110 ISW=1 C 120 I=NOB-2 C 130 DO 140 J=1,NA C B(NOB0+IS,J) C 140 B(I,J)=B(I,J) + B(NOB,J) C NOB=NOB0+IS C ISW=0 C GO TO 120 277 CONTINUE DO 2770 I=1,MX DO 2770 J=1,MX 2770 UM(I,J)=0.0D0 READ(15,*)IUM IF(IUM.NE.0)THEN DO 2772 I=1,NA-6 READ(15,*)J SUM=0.0D0 DO 2774 IK=1,J READ(15,*)IJ,UM(I,IJ) LV(IK)=IJ 2774 SUM=SUM+UM(I,IJ)*UM(I,IJ) DO 2773 IK=1,J 2773 UM(I,LV(IK))=UM(I,LV(IK))/sqrt(SUM) 2772 CONTINUE WRITE(*,*)' USER-DEFINED U MATRIX' M=NA-6 ELSE DO 2771 I=1,MX 2771 UM(I,I)=1.0D0 M=NOB0 ENDIF CLOSE(15) IF(IUM.NE.0)THEN DO 1556 I=1,NA-6 DO 1556 J=1,NA SUM=0.0D0 DO 1557 K=1,NOB0 1557 SUM=SUM+UM(I,K)*B(K,J) 1556 B0(I,J)=SUM ELSE DO 1558 I=1,NOB0 DO 1558 J=1,NA 1558 B0(I,J)=B(I,J) ENDIF DO 155 I=1,NOB0 DO 155 J=1,NA 155 IF(DABS(B(I,J)).LT.tol) B(I,J) = 0.0D0 WRITE(*,*)' OUTPUT FILE: B.MAT' OPEN(9,FILE='B.MAT') WRITE(9,901)M 901 format(I6,' Internal coordinates') WRITE(9,903)NA 903 format(I6,' Cartesian coordinates') WRITE(9,902) 902 format(' Transformation matrix, non-zero elements:') do 13 I=1,M WRITE(9,9011)I 9011 format(I6,' internal:') nn=0 do 14 J=1,NA 14 if(DABS(B0(I,J)).gt.tol)nn=nn+1 write(9,9012)nn 9012 format(I6) nn=0 do 13 J=1,NA if(DABS(B0(I,J)).gt.tol)then nn=nn+1 ia=(J+2)/3 ix=J-3*(ia-1) WRITE(9,9100)ia,ix,B0(I,J) 9100 format(2I5,F15.8) endif 13 continue CLOSE(9) CALL CTI WRITE(*,*) WRITE(*,*)' >> PROGRAM TERMINATED <<' STOP 190 IF(IER.EQ.1)WRITE(*,1130) 200 IF(ISCAN.EQ.0)WRITE(16,1170) ISCAN=ISCAN+1 IER=0 1070 FORMAT(/,1X,'coordinate ',I5,' of ',I5,':',/, 1 1X,'ILLEGAL CODE',I5,' CHOSEN') 1130 FORMAT(' ILLEGAL SPECIFICATION OF I, J, K, L, IX OR JX') 1170 FORMAT(/,1X,'*** PROGRAM WILL SCAN FOR FURTHER ERRORS IN DATA', $ ' DECK',/) END SUBROUTINE BOST C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR A BOND STRETCH C AS DEFINED BY WILSON. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=435) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,MX/3),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RIJ(3),RJK(3),RJL(3),EIJ(3),EJK(3),EKL(3) C IF(I.LE.0.OR.I.GT.NOAT)GO TO 30 IF(J.LE.0.OR.J.GT.NOAT)GO TO 30 IF(K.NE.0)GO TO 30 IF(L.NE.0)GO TO 30 IF(IX.NE.0)GO TO 30 IF(JX.NE.0)GO TO 30 DIJSQ=0.D0 DO 10 M=1,3 T=X(M,J)-X(M,I) RIJ(M)=T 10 DIJSQ=DIJSQ+T*T DIJ=DSQRT(DIJSQ) II=3*(I-1) JJ=3*(J-1) DO 20 M=1,3 T=RIJ(M) IF(DABS(T).LT.1.D-8)GO TO 20 T=T/DIJ B(NOB,II+M)=-T B(NOB,JJ+M)=T 20 CONTINUE RETURN 30 IER=1 WRITE(*,200) 200 FORMAT('BOST ERROR') RETURN END SUBROUTINE BEND C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS OF A VALENCE C ANGLE BENDING COORDINATE AS DEFINED BY WILSON. C I AND K ARE THE NUMBERS OF THE END ATOMS. C J IS THE NUMBER OF THE CENTRAL ATOM. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=435) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,MX/3),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RJI(3),RJK(3),RJL(3),EJI(3),EJK(3),EKL(3) C IF(I.LE.0.OR.I.GT.NOAT)GO TO 50 IF(J.LE.0.OR.J.GT.NOAT)GO TO 50 IF(K.LE.0.OR.K.GT.NOAT)GO TO 50 IF(L.NE.0)GO TO 50 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 50 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 50 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJISQ=0.D0 DJKSQ=0.D0 DXSQ=0.D0 DO 20 M=1,3 TP=X(M,J) T=X(M,I)-TP RJI(M)=T DJISQ=DJISQ+T*T T=X(M,K)-TP RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,JX)-X(M,IX) 20 DXSQ=DXSQ+T*T DJI=DSQRT(DJISQ) DJK=DSQRT(DJKSQ) DX=DSQRT(DXSQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 DOTJ=0.D0 DO 30 M=1,3 T=RJI(M)/DJI EJI(M)=T TP=RJK(M)/DJK EJK(M)=TP 30 DOTJ=DOTJ+T*TP IF(DABS(DOTJ).GT.0.99995D0)GO TO 60 SINJ=DSQRT(1.D0-DOTJ*DOTJ) II=3*(I-1) JJ=3*(J-1) KK=3*(K-1) DO 40 M=1,3 SMI=DX*(DOTJ*EJI(M)-EJK(M))/(DJI*SINJ) IF(DABS(SMI).GE.1.D-8)B(NOB,II+M)=SMI SMK=DX*(DOTJ*EJK(M)-EJI(M))/(DJK*SINJ) IF(DABS(SMK).GE.1.D-8)B(NOB,KK+M)=SMK SUM=SMI+SMK 40 IF(DABS(SUM).GE.1.D-8)B(NOB,JJ+M)=-SUM RETURN 50 IER=1 *C*C* ADDED BY KATH SPRING 1991 TO FIND INPUT ERROR * NEXT TWO LINES WRITE(*,200) 200 FORMAT('BEND ERROR') RETURN 60 IER=-1 WRITE(*,1000) RETURN C 1000 FORMAT(' I-J-K IS COLINEAR - USE LINEAR BEND') C END SUBROUTINE OPLA C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR AN OUT OF C PLANE WAGGING COORDINATE AS DEFINED BY DECIUS, MCINTOSH, MICHAELIAN C AND PETERSON. SUBROUTINE CODED BY M PETERSON, UNIV OF TORONTO. C I IS THE END ATOM (ATOM WAGGED WITH RESPECT TO J-K-L PLANE). C J IS THE APEX ATOM (ATOMS I, K AND L ARE ATTACHED TO J). C K AND L ARE THE ANCHOR ATOMS (DEFINE THE J-K-L PLANE). C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=435) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,MX/3),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RJI(3),RJK(3),RJL(3),EJI(3),EJK(3),EJL(3) DIMENSION C1(3) C IF(I.LE.0.OR.I.GT.NOAT)GO TO 60 IF(J.LE.0.OR.J.GT.NOAT)GO TO 60 IF(K.LE.0.OR.K.GT.NOAT)GO TO 60 IF(L.LE.0.OR.L.GT.NOAT)GO TO 60 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 60 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 60 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJISQ=0.D0 DJKSQ=0.D0 DJLSQ=0.D0 DXSQ=0.D0 DO 20 M=1,3 TP=X(M,J) T=X(M,I)-TP RJI(M)=T DJISQ=DJISQ+T*T T=X(M,K)-TP RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,L)-X(M,J) RJL(M)=T DJLSQ=DJLSQ+T*T T=X(M,JX)-X(M,IX) 20 DXSQ=DXSQ+T*T DJI=DSQRT(DJISQ) DJK=DSQRT(DJKSQ) DJL=DSQRT(DJLSQ) DX=DSQRT(DXSQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 COSI=0.D0 COSK=0.D0 COSL=0.D0 DO 30 M=1,3 T=RJI(M)/DJI EJI(M)=T TP=RJK(M)/DJK EJK(M)=TP TPP=RJL(M)/DJL EJL(M)=TPP COSI=COSI+TP*TPP COSK=COSK+T*TPP 30 COSL=COSL+T*TP IF(DABS(COSI).GT.0.99995D0)GO TO 70 SINSIN=1.D0-COSI*COSI SINI=DSQRT(SINSIN) C1(1)=EJK(2)*EJL(3)-EJK(3)*EJL(2) C1(2)=EJK(3)*EJL(1)-EJK(1)*EJL(3) C1(3)=EJK(1)*EJL(2)-EJK(2)*EJL(1) DOT=EJI(1)*C1(1)+EJI(2)*C1(2)+EJI(3)*C1(3) SINT=DOT/SINI IF(DABS(SINT).GT.0.00005D0)WRITE(*,1020) IF(DABS(SINT).GT.0.99995D0)GO TO 80 COST=DSQRT(1.D0-SINT*SINT) TANT=SINT/COST II=3*(I-1) JJ=3*(J-1) KK=3*(K-1) LL=3*(L-1) COSSIN=COST*SINI DO 50 M=1,3 T=C1(M)/COSSIN SMI=(T-TANT*EJI(M))/DJI IF(DABS(SMI).GE.1.D-8)B(NOB,II+M)=DX*SMI SMK=T*(COSI*COSK-COSL)/(SINSIN*DJK) IF(DABS(SMK).GE.1.D-8)B(NOB,KK+M)=DX*SMK SML=T*(COSI*COSL-COSK)/(SINSIN*DJL) IF(DABS(SML).GE.1.D-8)B(NOB,LL+M)=DX*SML SUM=SMI+SMK+SML 50 IF(DABS(SUM).GE.1.D-8)B(NOB,JJ+M)=-DX*SUM RETURN 60 IER=1 *C*C ADDED BY KATH, NEXT TWO LINES, SPRING 1991, LOOKING FOR ERRORS * WRITE(*,200) 200 FORMAT('OPLA ERROR') RETURN 70 IER=-1 WRITE(*,1000) RETURN 80 IER=-1 WRITE(*,1010) RETURN C 1000 FORMAT(' K-J-L IS COLINEAR (NO PLANE DEFINED FOR WAG OF I)') 1010 FORMAT(' I IS PERPENDICULAR TO J-K-L PLANE - USE VALENCE ANGLE 1 BENDS') 1020 FORMAT('+',86X,'*** WARNING: WAG OF A NON-PLANAR SYSTEM ***') C END SUBROUTINE TORS C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR TORSION AS DEFIN C BY R L HILDERBRANDT IN J MOLEC SPEC, 44, 599 (1972). C SUBROUTINE CODED BY M PETERSON, DEPT OF CHEMISTRY, UNIV OF TORONTO. C J AND K DEFINE THE BOND UNDER TORSION. C NI AND NL ARE THE NUMBER OF ATOMS ATTACHED TO J AND K RESPECTIVELY C (NI, NL <= 5). 2 DATA CARDS ARE READ: (1) CONTAINS NI ATOM NUMBERS C FOR THE I-TYPE ATOMS, AND (2) CONTAINS NL ATOM NUMBERS FOR THE L-TY C ATOMS (BOTH CARDS ARE IN 5I4 FORMAT). C C IATOM, LATOM: ATOM NUMBERS FOR THE I- AND L-TYPE ATOMS (SIZE: 5) C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=435) COMMON/BMAT/IC,NI,J,K,NL,IX,JX,NOAT,NOB,IER,X(3,MX/3),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RIJ(3),RJK(3),RLK(3),EIJ(3),EJK(3),ELK(3) DIMENSION CR(3),IATOM(5),LATOM(5),SJ(3),SK(3) C READ(15,*)(IATOM(I),I=1,NI) READ(15,*)(LATOM(L),L=1,NL) IF(NI.LE.0.OR.NI.GT.5)GO TO 110 IF(J.LE.0.OR.J.GT.NOAT)GO TO 110 IF(K.LE.0.OR.K.GT.NOAT)GO TO 110 IF(NL.LE.0.OR.NL.GT.5)GO TO 110 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 110 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 110 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJKSQ=0.D0 DXSQ=0.D0 DO 20 M=1,3 SJ(M)=0.D0 SK(M)=0.D0 T=X(M,K)-X(M,J) RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,JX)-X(M,IX) 20 DXSQ=DXSQ+T*T DJK=1.D0/DSQRT(DJKSQ) DX=DSQRT(DXSQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 DO 30 M=1,3 30 EJK(M)=RJK(M)*DJK JJ=3*(J-1) KK=3*(K-1) C C LOOP OVER THE I-TYPE ATOMS C DO 60 N=1,NI I=IATOM(N) IF(I.LE.0.OR.I.GT.NOAT)GO TO 110 DIJSQ=0.D0 DO 40 M=1,3 T=X(M,J)-X(M,I) RIJ(M)=T 40 DIJSQ=DIJSQ+T*T DIJ=1.D0/DSQRT(DIJSQ) COSJ=0.D0 DO 50 M=1,3 T=RIJ(M)*DIJ EIJ(M)=T 50 COSJ=COSJ-T*EJK(M) IF(DABS(COSJ).GT.0.99995D0)GO TO 120 SIN2J=(1.D0-COSJ*COSJ)*DFLOAT(NI) II=3*(I-1) CR(1)=EIJ(2)*EJK(3)-EIJ(3)*EJK(2) CR(2)=EIJ(3)*EJK(1)-EIJ(1)*EJK(3) CR(3)=EIJ(1)*EJK(2)-EIJ(2)*EJK(1) DO 60 M=1,3 T=CR(M)/SIN2J SMI=T*DIJ IF(DABS(SMI).GE.1.D-8)B(NOB,II+M)=-DX*SMI SMK=T*COSJ*DJK SK(M)=SK(M)+SMK SMJ=SMI-SMK 60 SJ(M)=SJ(M)+SMJ C C LOOP OVER THE L-TYPE ATOMS C DO 90 N=1,NL L=LATOM(N) IF(L.LE.0.OR.L.GT.NOAT)GO TO 110 DLKSQ=0.D0 DO 70 M=1,3 T=X(M,K)-X(M,L) RLK(M)=T 70 DLKSQ=DLKSQ+T*T DLK=1.D0/DSQRT(DLKSQ) COSK=0.D0 DO 80 M=1,3 T=RLK(M)*DLK ELK(M)=T 80 COSK=COSK+EJK(M)*T IF(DABS(COSK).GT.0.99995D0)GO TO 120 SIN2K=(1.D0-COSK*COSK)*DFLOAT(NL) LL=3*(L-1) CR(1)=ELK(3)*EJK(2)-ELK(2)*EJK(3) CR(2)=ELK(1)*EJK(3)-ELK(3)*EJK(1) CR(3)=ELK(2)*EJK(1)-ELK(1)*EJK(2) DO 90 M=1,3 T=CR(M)/SIN2K SML=T*DLK IF(DABS(SML).GE.1.D-8)B(NOB,LL+M)=-DX*SML SMJ=T*COSK*DJK SJ(M)=SJ(M)+SMJ SMK=SML-SMJ 90 SK(M)=SK(M)+SMK DO 100 M=1,3 SMJ=SJ(M) IF(DABS(SMJ).GE.1.D-8)B(NOB,JJ+M)=SMJ*DX SMK=SK(M) 100 IF(DABS(SMK).GE.1.D-8)B(NOB,KK+M)=SMK*DX RETURN 110 IER=1 WRITE(*,200) I,L 200 FORMAT('TORS ERROR AT I=',I2,' AND L=',I2) RETURN 120 IER=-1 WRITE(*,1030) RETURN 1030 FORMAT(' I-J-K OR J-K-L IS COLINEAR (NO TORSION POSSIBLE)') END SUBROUTINE LIBE C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR A LINEAR BEND C OR FOR A PAIR OF PERPENDICULAR LINEAR BENDS. C I AND K ARE THE END ATOMS. C J IS THE CENTRAL ATOM. C C A GIVES THE CARTESIAN COORDINATES OF A POINT IN SPACE, SUCH C THAT THE VECTOR FROM ATOM J TO POINT A IS PERPENDICULAR TO C THE LINE I-J-K AND SERVES TO ORIENT THE COORDINATES IN SPACE. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=435) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,MX/3),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RJI(3),RJK(3),EJK(3),UP(3),UN(3),UNIT(3) DIMENSION A(3) C READ(15,1000)A IF(I.LE.0.OR.I.GT.NOAT)GO TO 60 IF(J.LE.0.OR.J.GT.NOAT)GO TO 60 IF(K.LE.0.OR.K.GT.NOAT)GO TO 60 IF(L.NE.0)GO TO 60 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 60 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 60 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJISQ=0.D0 DJKSQ=0.D0 DXSQ=0.D0 DJASQ=0.D0 DO 20 M=1,3 TP=X(M,J) T=X(M,I)-TP RJI(M)=T DJISQ=DJISQ+T*T T=X(M,K)-TP RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,JX)-X(M,IX) DXSQ=DXSQ+T*T T=A(M)-TP UN(M)=T 20 DJASQ=DJASQ+T*T DJI=DSQRT(DJISQ) DJK=DSQRT(DJKSQ) DX=DSQRT(DXSQ) DJA=DSQRT(DJASQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 DOTJ=0.D0 DOTP=0.D0 DO 30 M=1,3 T=RJI(M)/DJI TP=RJK(M)/DJK EJK(M)=TP DOTJ=DOTJ+T*TP TP=UN(M)/DJA UNIT(M)=TP 30 DOTP=DOTP+T*TP TEST=DABS(DOTJ)-1.D0 IF(DABS(TEST).GT.0.00005D0)GO TO 70 IF(DABS(DOTP).GT.0.00005D0)GO TO 80 II=3*(I-1) JJ=3*(J-1) KK=3*(K-1) DO 40 M=1,3 T=UNIT(M) IF(DABS(T).LT.1.D-8)GO TO 40 T=-DX*T SMI=T/DJI B(NOB,II+M)=SMI SMK=T/DJK B(NOB,KK+M)=SMK B(NOB,JJ+M)=-SMI-SMK 40 CONTINUE IF(IC.EQ.6)RETURN NOB=NOB0+IS UP(1)=EJK(2)*UNIT(3)-EJK(3)*UNIT(2) UP(2)=EJK(3)*UNIT(1)-EJK(1)*UNIT(3) UP(3)=EJK(1)*UNIT(2)-EJK(2)*UNIT(1) DO 50 M=1,3 T=UP(M) IF(DABS(T).LT.1.D-8)GO TO 50 T=-DX*T SMI=T/DJI B(NOB,II+M)=SMI SMK=T/DJK B(NOB,KK+M)=SMK B(NOB,JJ+M)=-SMI-SMK 50 CONTINUE RETURN 60 IER=1 RETURN 70 IER=-1 WRITE(*,1020) RETURN 80 IER=-1 WRITE(*,1030) RETURN 1000 FORMAT(3D23.16) 1020 FORMAT(' I-J-K NOT COLINEAR - USE VALENCE ANGLE BEND') 1030 FORMAT(' ATOM A NOT PERPENDICULAR TO I-J-K AT J') END SUBROUTINE CTI IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (MX=435) DIMENSION FINT(MX,MX),B(MX,MX),FCAR(MX,MX) allocatable::tem(:),am(:,:),BV(:,:),w(:,:) real*8,allocatable::WORK(:) integer*4,allocatable::IPIV(:) integer*4 time WRITE(*,*)' INTRISIC FORCE FIELD GENERATION' open(9,file='B.MAT') read(9,*)nob0 read(9,*)na do 14 i=1,nob0 do 14 j=1,na 14 b(i,j)=0.0d0 read(9,*) do 13 i=1,nob0 read(9,*) read(9,*)nn do 13 ii=1,nn 13 read(9,*)ia,ix,b(i,ix+3*(ia-1)) close(9) write(6,*)'B.MAT read, ',na,' atomic ',nob0,' internal' NAT3=NA N=NAT3 call readff(MX,N,FCAR) allocate(am(nob0*nob0,nob0*nob0),tem(nob0*nob0)) alpha=1.0d-6 c c make AM = BBBB do 1 i1=1,nob0 do 1 i2=1,nob0 tem(i1+(i2-1)*nob0)=0.0d0 do 1 lc=1,na 1 tem(i1+(i2-1)*nob0)=tem(i1+(i2-1)*nob0)+b(i1,lc)*b(i2,lc) write(6,*)'t done' c do 5 i1=1,nob0 do 5 i2=1,nob0 i=i1+(i2-1)*nob0 do 5 j1=1,nob0 do 5 j2=1,nob0 j=j1+(j2-1)*nob0 am(i,j)=tem(i1+(j1-1)*nob0)*tem(i2+(j2-1)*nob0) 5 if(i.eq.j)am(i,j)=am(i,j)+alpha write(6,*)'a done' deallocate(tem) it0=time() NT=nob0**2 write(6,607)NT,NT 607 format(' A X =B, dimensions ',i9,' x ',i9) c Form BV = BBF: allocate(w(nob0,N),BV(NT,1)) c first index do 3 i1=1,nob0 do 3 ic=1,N w(i1,ic)=0.0d0 do 3 lc=1,N 3 w(i1,ic)=w(i1,ic)+FCAR(ic,lc)*b(i1,lc) c second index do 7 i1=1,nob0 do 7 i2=1,nob0 BV(i1+(i2-1)*nob0,1)=0.0d0 do 7 lc=1,N 7 BV(i1+(i2-1)*nob0,1)=BV(i1+(i2-1)*nob0,1)+w(i1,lc)*b(i2,lc) write(6,*)'c done' deallocate(w) c linear solver without matrix inversion: A FI = BV INFO=0 allocate(WORK(1),IPIV(NT)) LWORK=-1 c find dimension of the working buffer: call DSYTRF('U', NT,am, NT, IPIV, WORK, LWORK, INFO ) LWORK=INT(WORK(1)) if(INFO.eq.0)write(6,*)'UL success',LWORK deallocate(WORK) allocate(WORK(LWORK)) c factorize the matrix: call DSYTRF('U', NT, am, NT, IPIV, WORK, LWORK, INFO ) INFO=0 c find the solution, put in BV: call DSYTRS('U',NT, 1, am, NT, IPIV, BV, NT, INFO ) if(INFO.eq.0)write(6,*)'success' write(6,609)time()-it0 609 format('inversion time ',i20) deallocate(am) c rewrite:BV to FINT do 4 i=1,nob0 do 4 j=1,nob0 4 FINT(i,j)=BV(i+(j-1)*nob0,1) WRITE(*,*)'INTY.FC : OUTPUT WITH THE INTRISIC FORCE CONSTANTS' WRITE(*,*)'INTY.SCR : BINARY' OPEN(2,FILE='INTY.FC') WRITE(2,*)' Force constants in internal coordinates:' CALL WRITEFF(MX,nob0,FINT) close(2) OPEN(29,FILE='INTY.SCR',FORM='UNFORMATTED') NINT=nob0 LTR=(NINT*(NINT+1))/2 WRITE(29) LTR,((FINT(I,J),I=1,J),J=1,nob0) CLOSE(29) write(*,*)'written' RETURN END SUBROUTINE WRITEFF(MX,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX,MX) N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N write(2,18)(ij,ij=N1,N3) 18 FORMAT(4x,5I14) DO 130 LN=N1,N 130 WRITE(2,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) RETURN END *> \brief \b DSYTRS * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download DSYTRS + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE DSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO * INTEGER INFO, LDA, LDB, N, NRHS * .. * .. Array Arguments .. * INTEGER IPIV( * ) * DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DSYTRS solves a system of linear equations A*X = B with a real *> symmetric matrix A using the factorization A = U*D*U**T or *> A = L*D*L**T computed by DSYTRF. *> \endverbatim * * Arguments: * ========== * *> \param[in] UPLO *> \verbatim *> UPLO is CHARACTER*1 *> Specifies whether the details of the factorization are stored *> as an upper or lower triangular matrix. *> = 'U': Upper triangular, form is A = U*D*U**T; *> = 'L': Lower triangular, form is A = L*D*L**T. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> The order of the matrix A. N >= 0. *> \endverbatim *> *> \param[in] NRHS *> \verbatim *> NRHS is INTEGER *> The number of right hand sides, i.e., the number of columns *> of the matrix B. NRHS >= 0. *> \endverbatim *> *> \param[in] A *> \verbatim *> A is DOUBLE PRECISION array, dimension (LDA,N) *> The block diagonal matrix D and the multipliers used to *> obtain the factor U or L as computed by DSYTRF. *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is INTEGER *> The leading dimension of the array A. LDA >= max(1,N). *> \endverbatim *> *> \param[in] IPIV *> \verbatim *> IPIV is INTEGER array, dimension (N) *> Details of the interchanges and the block structure of D *> as determined by DSYTRF. *> \endverbatim *> *> \param[in,out] B *> \verbatim *> B is DOUBLE PRECISION array, dimension (LDB,NRHS) *> On entry, the right hand side matrix B. *> On exit, the solution matrix X. *> \endverbatim *> *> \param[in] LDB *> \verbatim *> LDB is INTEGER *> The leading dimension of the array B. LDB >= max(1,N). *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2011 * *> \ingroup doubleSYcomputational * * ===================================================================== SUBROUTINE DSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) * * -- LAPACK computational routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, LDA, LDB, N, NRHS * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL UPPER INTEGER J, K, KP DOUBLE PRECISION AK, AKM1, AKM1K, BK, BKM1, DENOM * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DGEMV, DGER, DSCAL, DSWAP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( NRHS.LT.0 ) THEN INFO = -3 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -5 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -8 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DSYTRS', -INFO ) RETURN END IF * * Quick return if possible * IF( N.EQ.0 .OR. NRHS.EQ.0 ) $ RETURN * IF( UPPER ) THEN * * Solve A*X = B, where A = U*D*U**T. * * First solve U*D*X = B, overwriting B with X. * * K is the main loop index, decreasing from N to 1 in steps of * 1 or 2, depending on the size of the diagonal blocks. * K = N 10 CONTINUE * * If K < 1, exit from loop. * IF( K.LT.1 ) $ GO TO 30 * IF( IPIV( K ).GT.0 ) THEN * * 1 x 1 diagonal block * * Interchange rows K and IPIV(K). * KP = IPIV( K ) IF( KP.NE.K ) $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) * * Multiply by inv(U(K)), where U(K) is the transformation * stored in column K of A. * CALL DGER( K-1, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB, $ B( 1, 1 ), LDB ) * * Multiply by the inverse of the diagonal block. * CALL DSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB ) K = K - 1 ELSE * * 2 x 2 diagonal block * * Interchange rows K-1 and -IPIV(K). * KP = -IPIV( K ) IF( KP.NE.K-1 ) $ CALL DSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ), LDB ) * * Multiply by inv(U(K)), where U(K) is the transformation * stored in columns K-1 and K of A. * CALL DGER( K-2, NRHS, -ONE, A( 1, K ), 1, B( K, 1 ), LDB, $ B( 1, 1 ), LDB ) CALL DGER( K-2, NRHS, -ONE, A( 1, K-1 ), 1, B( K-1, 1 ), $ LDB, B( 1, 1 ), LDB ) * * Multiply by the inverse of the diagonal block. * AKM1K = A( K-1, K ) AKM1 = A( K-1, K-1 ) / AKM1K AK = A( K, K ) / AKM1K DENOM = AKM1*AK - ONE DO 20 J = 1, NRHS BKM1 = B( K-1, J ) / AKM1K BK = B( K, J ) / AKM1K B( K-1, J ) = ( AK*BKM1-BK ) / DENOM B( K, J ) = ( AKM1*BK-BKM1 ) / DENOM 20 CONTINUE K = K - 2 END IF * GO TO 10 30 CONTINUE * * Next solve U**T *X = B, overwriting B with X. * * K is the main loop index, increasing from 1 to N in steps of * 1 or 2, depending on the size of the diagonal blocks. * K = 1 40 CONTINUE * * If K > N, exit from loop. * IF( K.GT.N ) $ GO TO 50 * IF( IPIV( K ).GT.0 ) THEN * * 1 x 1 diagonal block * * Multiply by inv(U**T(K)), where U(K) is the transformation * stored in column K of A. * CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, A( 1, K ), $ 1, ONE, B( K, 1 ), LDB ) * * Interchange rows K and IPIV(K). * KP = IPIV( K ) IF( KP.NE.K ) $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) K = K + 1 ELSE * * 2 x 2 diagonal block * * Multiply by inv(U**T(K+1)), where U(K+1) is the transformation * stored in columns K and K+1 of A. * CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, A( 1, K ), $ 1, ONE, B( K, 1 ), LDB ) CALL DGEMV( 'Transpose', K-1, NRHS, -ONE, B, LDB, $ A( 1, K+1 ), 1, ONE, B( K+1, 1 ), LDB ) * * Interchange rows K and -IPIV(K). * KP = -IPIV( K ) IF( KP.NE.K ) $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) K = K + 2 END IF * GO TO 40 50 CONTINUE * ELSE * * Solve A*X = B, where A = L*D*L**T. * * First solve L*D*X = B, overwriting B with X. * * K is the main loop index, increasing from 1 to N in steps of * 1 or 2, depending on the size of the diagonal blocks. * K = 1 60 CONTINUE * * If K > N, exit from loop. * IF( K.GT.N ) $ GO TO 80 * IF( IPIV( K ).GT.0 ) THEN * * 1 x 1 diagonal block * * Interchange rows K and IPIV(K). * KP = IPIV( K ) IF( KP.NE.K ) $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) * * Multiply by inv(L(K)), where L(K) is the transformation * stored in column K of A. * IF( K.LT.N ) $ CALL DGER( N-K, NRHS, -ONE, A( K+1, K ), 1, B( K, 1 ), $ LDB, B( K+1, 1 ), LDB ) * * Multiply by the inverse of the diagonal block. * CALL DSCAL( NRHS, ONE / A( K, K ), B( K, 1 ), LDB ) K = K + 1 ELSE * * 2 x 2 diagonal block * * Interchange rows K+1 and -IPIV(K). * KP = -IPIV( K ) IF( KP.NE.K+1 ) $ CALL DSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ), LDB ) * * Multiply by inv(L(K)), where L(K) is the transformation * stored in columns K and K+1 of A. * IF( K.LT.N-1 ) THEN CALL DGER( N-K-1, NRHS, -ONE, A( K+2, K ), 1, B( K, 1 ), $ LDB, B( K+2, 1 ), LDB ) CALL DGER( N-K-1, NRHS, -ONE, A( K+2, K+1 ), 1, $ B( K+1, 1 ), LDB, B( K+2, 1 ), LDB ) END IF * * Multiply by the inverse of the diagonal block. * AKM1K = A( K+1, K ) AKM1 = A( K, K ) / AKM1K AK = A( K+1, K+1 ) / AKM1K DENOM = AKM1*AK - ONE DO 70 J = 1, NRHS BKM1 = B( K, J ) / AKM1K BK = B( K+1, J ) / AKM1K B( K, J ) = ( AK*BKM1-BK ) / DENOM B( K+1, J ) = ( AKM1*BK-BKM1 ) / DENOM 70 CONTINUE K = K + 2 END IF * GO TO 60 80 CONTINUE * * Next solve L**T *X = B, overwriting B with X. * * K is the main loop index, decreasing from N to 1 in steps of * 1 or 2, depending on the size of the diagonal blocks. * K = N 90 CONTINUE * * If K < 1, exit from loop. * IF( K.LT.1 ) $ GO TO 100 * IF( IPIV( K ).GT.0 ) THEN * * 1 x 1 diagonal block * * Multiply by inv(L**T(K)), where L(K) is the transformation * stored in column K of A. * IF( K.LT.N ) $ CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ), $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB ) * * Interchange rows K and IPIV(K). * KP = IPIV( K ) IF( KP.NE.K ) $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) K = K - 1 ELSE * * 2 x 2 diagonal block * * Multiply by inv(L**T(K-1)), where L(K-1) is the transformation * stored in columns K-1 and K of A. * IF( K.LT.N ) THEN CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ), $ LDB, A( K+1, K ), 1, ONE, B( K, 1 ), LDB ) CALL DGEMV( 'Transpose', N-K, NRHS, -ONE, B( K+1, 1 ), $ LDB, A( K+1, K-1 ), 1, ONE, B( K-1, 1 ), $ LDB ) END IF * * Interchange rows K and -IPIV(K). * KP = -IPIV( K ) IF( KP.NE.K ) $ CALL DSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) K = K - 2 END IF * GO TO 90 100 CONTINUE END IF * RETURN * * End of DSYTRS * END *> \brief \b DSYTRF * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download DSYTRF + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE DSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO * INTEGER INFO, LDA, LWORK, N * .. * .. Array Arguments .. * INTEGER IPIV( * ) * DOUBLE PRECISION A( LDA, * ), WORK( * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DSYTRF computes the factorization of a real symmetric matrix A using *> the Bunch-Kaufman diagonal pivoting method. The form of the *> factorization is *> *> A = U*D*U**T or A = L*D*L**T *> *> where U (or L) is a product of permutation and unit upper (lower) *> triangular matrices, and D is symmetric and block diagonal with *> 1-by-1 and 2-by-2 diagonal blocks. *> *> This is the blocked version of the algorithm, calling Level 3 BLAS. *> \endverbatim * * Arguments: * ========== * *> \param[in] UPLO *> \verbatim *> UPLO is CHARACTER*1 *> = 'U': Upper triangle of A is stored; *> = 'L': Lower triangle of A is stored. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> The order of the matrix A. N >= 0. *> \endverbatim *> *> \param[in,out] A *> \verbatim *> A is DOUBLE PRECISION array, dimension (LDA,N) *> On entry, the symmetric matrix A. If UPLO = 'U', the leading *> N-by-N upper triangular part of A contains the upper *> triangular part of the matrix A, and the strictly lower *> triangular part of A is not referenced. If UPLO = 'L', the *> leading N-by-N lower triangular part of A contains the lower *> triangular part of the matrix A, and the strictly upper *> triangular part of A is not referenced. *> *> On exit, the block diagonal matrix D and the multipliers used *> to obtain the factor U or L (see below for further details). *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is INTEGER *> The leading dimension of the array A. LDA >= max(1,N). *> \endverbatim *> *> \param[out] IPIV *> \verbatim *> IPIV is INTEGER array, dimension (N) *> Details of the interchanges and the block structure of D. *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were *> interchanged and D(k,k) is a 1-by-1 diagonal block. *> If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and *> columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) *> is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) = *> IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were *> interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block. *> \endverbatim *> *> \param[out] WORK *> \verbatim *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. *> \endverbatim *> *> \param[in] LWORK *> \verbatim *> LWORK is INTEGER *> The length of WORK. LWORK >=1. For best performance *> LWORK >= N*NB, where NB is the block size 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. *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> > 0: if INFO = i, D(i,i) is exactly zero. The factorization *> has been completed, but the block diagonal matrix D is *> exactly singular, and division by zero will occur if it *> is used to solve a system of equations. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2011 * *> \ingroup doubleSYcomputational * *> \par Further Details: * ===================== *> *> \verbatim *> *> If UPLO = 'U', then A = U*D*U**T, where *> U = P(n)*U(n)* ... *P(k)U(k)* ..., *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such *> that if the diagonal block D(k) is of order s (s = 1 or 2), then *> *> ( I v 0 ) k-s *> U(k) = ( 0 I 0 ) s *> ( 0 0 I ) n-k *> k-s s n-k *> *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), *> and A(k,k), and v overwrites A(1:k-2,k-1:k). *> *> If UPLO = 'L', then A = L*D*L**T, where *> L = P(1)*L(1)* ... *P(k)*L(k)* ..., *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such *> that if the diagonal block D(k) is of order s (s = 1 or 2), then *> *> ( I 0 0 ) k-1 *> L(k) = ( 0 I 0 ) s *> ( 0 v I ) n-k-s+1 *> k-1 s n-k-s+1 *> *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). *> \endverbatim *> * ===================================================================== SUBROUTINE DSYTRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO ) * * -- LAPACK computational routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, LDA, LWORK, N * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ), WORK( * ) * .. * * ===================================================================== * * .. Local Scalars .. LOGICAL LQUERY, UPPER INTEGER IINFO, IWS, J, K, KB, LDWORK, LWKOPT, NB, NBMIN * .. * .. External Functions .. LOGICAL LSAME INTEGER ILAENV EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. EXTERNAL DLASYF, DSYTF2, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) LQUERY = ( LWORK.EQ.-1 ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN INFO = -7 END IF * IF( INFO.EQ.0 ) THEN * * Determine the block size * NB = ILAENV( 1, 'DSYTRF', UPLO, N, -1, -1, -1 ) LWKOPT = N*NB WORK( 1 ) = LWKOPT END IF * IF( INFO.NE.0 ) THEN CALL XERBLA( 'DSYTRF', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF * NBMIN = 2 LDWORK = N IF( NB.GT.1 .AND. NB.LT.N ) THEN IWS = LDWORK*NB IF( LWORK.LT.IWS ) THEN NB = MAX( LWORK / LDWORK, 1 ) NBMIN = MAX( 2, ILAENV( 2, 'DSYTRF', UPLO, N, -1, -1, -1 ) ) END IF ELSE IWS = 1 END IF IF( NB.LT.NBMIN ) $ NB = N * IF( UPPER ) THEN * * Factorize A as U*D*U**T using the upper triangle of A * * K is the main loop index, decreasing from N to 1 in steps of * KB, where KB is the number of columns factorized by DLASYF; * KB is either NB or NB-1, or K for the last block * K = N 10 CONTINUE * * If K < 1, exit from loop * IF( K.LT.1 ) $ GO TO 40 * IF( K.GT.NB ) THEN * * Factorize columns k-kb+1:k of A and use blocked code to * update columns 1:k-kb * CALL DLASYF( UPLO, K, NB, KB, A, LDA, IPIV, WORK, LDWORK, $ IINFO ) ELSE * * Use unblocked code to factorize columns 1:k of A * CALL DSYTF2( UPLO, K, A, LDA, IPIV, IINFO ) KB = K END IF * * Set INFO on the first occurrence of a zero pivot * IF( INFO.EQ.0 .AND. IINFO.GT.0 ) $ INFO = IINFO * * Decrease K and return to the start of the main loop * K = K - KB GO TO 10 * ELSE * * Factorize A as L*D*L**T using the lower triangle of A * * K is the main loop index, increasing from 1 to N in steps of * KB, where KB is the number of columns factorized by DLASYF; * KB is either NB or NB-1, or N-K+1 for the last block * K = 1 20 CONTINUE * * If K > N, exit from loop * IF( K.GT.N ) $ GO TO 40 * IF( K.LE.N-NB ) THEN * * Factorize columns k:k+kb-1 of A and use blocked code to * update columns k+kb:n * CALL DLASYF( UPLO, N-K+1, NB, KB, A( K, K ), LDA, IPIV( K ), $ WORK, LDWORK, IINFO ) ELSE * * Use unblocked code to factorize columns k:n of A * CALL DSYTF2( UPLO, N-K+1, A( K, K ), LDA, IPIV( K ), IINFO ) KB = N - K + 1 END IF * * Set INFO on the first occurrence of a zero pivot * IF( INFO.EQ.0 .AND. IINFO.GT.0 ) $ INFO = IINFO + K - 1 * * Adjust IPIV * DO 30 J = K, K + KB - 1 IF( IPIV( J ).GT.0 ) THEN IPIV( J ) = IPIV( J ) + K - 1 ELSE IPIV( J ) = IPIV( J ) - K + 1 END IF 30 CONTINUE * * Increase K and return to the start of the main loop * K = K + KB GO TO 20 * END IF * 40 CONTINUE WORK( 1 ) = LWKOPT RETURN * * End of DSYTRF * END *> \brief \b DLASYF computes a partial factorization of a real symmetric matrix using the Bunch-Kaufman diagonal pivoting method. * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download DLASYF + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO * INTEGER INFO, KB, LDA, LDW, N, NB * .. * .. Array Arguments .. * INTEGER IPIV( * ) * DOUBLE PRECISION A( LDA, * ), W( LDW, * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DLASYF computes a partial factorization of a real symmetric matrix A *> using the Bunch-Kaufman diagonal pivoting method. The partial *> factorization has the form: *> *> A = ( I U12 ) ( A11 0 ) ( I 0 ) if UPLO = 'U', or: *> ( 0 U22 ) ( 0 D ) ( U12**T U22**T ) *> *> A = ( L11 0 ) ( D 0 ) ( L11**T L21**T ) if UPLO = 'L' *> ( L21 I ) ( 0 A22 ) ( 0 I ) *> *> where the order of D is at most NB. The actual order is returned in *> the argument KB, and is either NB or NB-1, or N if N <= NB. *> *> DLASYF is an auxiliary routine called by DSYTRF. It uses blocked code *> (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or *> A22 (if UPLO = 'L'). *> \endverbatim * * Arguments: * ========== * *> \param[in] UPLO *> \verbatim *> UPLO is CHARACTER*1 *> Specifies whether the upper or lower triangular part of the *> symmetric matrix A is stored: *> = 'U': Upper triangular *> = 'L': Lower triangular *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> The order of the matrix A. N >= 0. *> \endverbatim *> *> \param[in] NB *> \verbatim *> NB is INTEGER *> The maximum number of columns of the matrix A that should be *> factored. NB should be at least 2 to allow for 2-by-2 pivot *> blocks. *> \endverbatim *> *> \param[out] KB *> \verbatim *> KB is INTEGER *> The number of columns of A that were actually factored. *> KB is either NB-1 or NB, or N if N <= NB. *> \endverbatim *> *> \param[in,out] A *> \verbatim *> A is DOUBLE PRECISION array, dimension (LDA,N) *> On entry, the symmetric matrix A. If UPLO = 'U', the leading *> n-by-n upper triangular part of A contains the upper *> triangular part of the matrix A, and the strictly lower *> triangular part of A is not referenced. If UPLO = 'L', the *> leading n-by-n lower triangular part of A contains the lower *> triangular part of the matrix A, and the strictly upper *> triangular part of A is not referenced. *> On exit, A contains details of the partial factorization. *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is INTEGER *> The leading dimension of the array A. LDA >= max(1,N). *> \endverbatim *> *> \param[out] IPIV *> \verbatim *> IPIV is INTEGER array, dimension (N) *> Details of the interchanges and the block structure of D. *> *> If UPLO = 'U': *> Only the last KB elements of IPIV are set. *> *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were *> interchanged and D(k,k) is a 1-by-1 diagonal block. *> *> If IPIV(k) = IPIV(k-1) < 0, then rows and columns *> k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) *> is a 2-by-2 diagonal block. *> *> If UPLO = 'L': *> Only the first KB elements of IPIV are set. *> *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were *> interchanged and D(k,k) is a 1-by-1 diagonal block. *> *> If IPIV(k) = IPIV(k+1) < 0, then rows and columns *> k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1) *> is a 2-by-2 diagonal block. *> \endverbatim *> *> \param[out] W *> \verbatim *> W is DOUBLE PRECISION array, dimension (LDW,NB) *> \endverbatim *> *> \param[in] LDW *> \verbatim *> LDW is INTEGER *> The leading dimension of the array W. LDW >= max(1,N). *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is INTEGER *> = 0: successful exit *> > 0: if INFO = k, D(k,k) is exactly zero. The factorization *> has been completed, but the block diagonal matrix D is *> exactly singular. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2013 * *> \ingroup doubleSYcomputational * *> \par Contributors: * ================== *> *> \verbatim *> *> November 2013, Igor Kozachenko, *> Computer Science Division, *> University of California, Berkeley *> \endverbatim * * ===================================================================== SUBROUTINE DLASYF( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, INFO ) * * -- LAPACK computational routine (version 3.5.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2013 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, KB, LDA, LDW, N, NB * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ), W( LDW, * ) * .. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) DOUBLE PRECISION EIGHT, SEVTEN PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 ) * .. * .. Local Scalars .. INTEGER IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP, $ KSTEP, KW DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D21, D22, R1, $ ROWMAX, T * .. * .. External Functions .. LOGICAL LSAME INTEGER IDAMAX EXTERNAL LSAME, IDAMAX * .. * .. External Subroutines .. EXTERNAL DCOPY, DGEMM, DGEMV, DSCAL, DSWAP * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. * .. Executable Statements .. * INFO = 0 * * Initialize ALPHA for use in choosing pivot block size. * ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT * IF( LSAME( UPLO, 'U' ) ) THEN * * Factorize the trailing columns of A using the upper triangle * of A and working backwards, and compute the matrix W = U12*D * for use in updating A11 * * K is the main loop index, decreasing from N in steps of 1 or 2 * * KW is the column of W which corresponds to column K of A * K = N 10 CONTINUE KW = NB + K - N * * Exit from loop * IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 ) $ GO TO 30 * * Copy column K of A to column KW of W and update it * CALL DCOPY( K, A( 1, K ), 1, W( 1, KW ), 1 ) IF( K.LT.N ) $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ), LDA, $ W( K, KW+1 ), LDW, ONE, W( 1, KW ), 1 ) * KSTEP = 1 * * Determine rows and columns to be interchanged and whether * a 1-by-1 or 2-by-2 pivot block will be used * ABSAKK = ABS( W( K, KW ) ) * * IMAX is the row-index of the largest off-diagonal element in * column K, and COLMAX is its absolute value. * Determine both COLMAX and IMAX. * IF( K.GT.1 ) THEN IMAX = IDAMAX( K-1, W( 1, KW ), 1 ) COLMAX = ABS( W( IMAX, KW ) ) ELSE COLMAX = ZERO END IF * IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN * * Column K is zero or underflow: set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K KP = K ELSE IF( ABSAKK.GE.ALPHA*COLMAX ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE * * Copy column IMAX to column KW-1 of W and update it * CALL DCOPY( IMAX, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 ) CALL DCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA, $ W( IMAX+1, KW-1 ), 1 ) IF( K.LT.N ) $ CALL DGEMV( 'No transpose', K, N-K, -ONE, A( 1, K+1 ), $ LDA, W( IMAX, KW+1 ), LDW, ONE, $ W( 1, KW-1 ), 1 ) * * JMAX is the column-index of the largest off-diagonal * element in row IMAX, and ROWMAX is its absolute value * JMAX = IMAX + IDAMAX( K-IMAX, W( IMAX+1, KW-1 ), 1 ) ROWMAX = ABS( W( JMAX, KW-1 ) ) IF( IMAX.GT.1 ) THEN JMAX = IDAMAX( IMAX-1, W( 1, KW-1 ), 1 ) ROWMAX = MAX( ROWMAX, ABS( W( JMAX, KW-1 ) ) ) END IF * IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE IF( ABS( W( IMAX, KW-1 ) ).GE.ALPHA*ROWMAX ) THEN * * interchange rows and columns K and IMAX, use 1-by-1 * pivot block * KP = IMAX * * copy column KW-1 of W to column KW of W * CALL DCOPY( K, W( 1, KW-1 ), 1, W( 1, KW ), 1 ) ELSE * * interchange rows and columns K-1 and IMAX, use 2-by-2 * pivot block * KP = IMAX KSTEP = 2 END IF END IF * * ============================================================ * * KK is the column of A where pivoting step stopped * KK = K - KSTEP + 1 * * KKW is the column of W which corresponds to column KK of A * KKW = NB + KK - N * * Interchange rows and columns KP and KK. * Updated column KP is already stored in column KKW of W. * IF( KP.NE.KK ) THEN * * Copy non-updated column KK to column KP of submatrix A * at step K. No need to copy element into column K * (or K and K-1 for 2-by-2 pivot) of A, since these columns * will be later overwritten. * A( KP, KP ) = A( KK, KK ) CALL DCOPY( KK-1-KP, A( KP+1, KK ), 1, A( KP, KP+1 ), $ LDA ) IF( KP.GT.1 ) $ CALL DCOPY( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 ) * * Interchange rows KK and KP in last K+1 to N columns of A * (columns K (or K and K-1 for 2-by-2 pivot) of A will be * later overwritten). Interchange rows KK and KP * in last KKW to NB columns of W. * IF( K.LT.N ) $ CALL DSWAP( N-K, A( KK, K+1 ), LDA, A( KP, K+1 ), $ LDA ) CALL DSWAP( N-KK+1, W( KK, KKW ), LDW, W( KP, KKW ), $ LDW ) END IF * IF( KSTEP.EQ.1 ) THEN * * 1-by-1 pivot block D(k): column kw of W now holds * * W(kw) = U(k)*D(k), * * where U(k) is the k-th column of U * * Store subdiag. elements of column U(k) * and 1-by-1 block D(k) in column k of A. * NOTE: Diagonal element U(k,k) is a UNIT element * and not stored. * A(k,k) := D(k,k) = W(k,kw) * A(1:k-1,k) := U(1:k-1,k) = W(1:k-1,kw)/D(k,k) * CALL DCOPY( K, W( 1, KW ), 1, A( 1, K ), 1 ) R1 = ONE / A( K, K ) CALL DSCAL( K-1, R1, A( 1, K ), 1 ) * ELSE * * 2-by-2 pivot block D(k): columns kw and kw-1 of W now hold * * ( W(kw-1) W(kw) ) = ( U(k-1) U(k) )*D(k) * * where U(k) and U(k-1) are the k-th and (k-1)-th columns * of U * * Store U(1:k-2,k-1) and U(1:k-2,k) and 2-by-2 * block D(k-1:k,k-1:k) in columns k-1 and k of A. * NOTE: 2-by-2 diagonal block U(k-1:k,k-1:k) is a UNIT * block and not stored. * A(k-1:k,k-1:k) := D(k-1:k,k-1:k) = W(k-1:k,kw-1:kw) * A(1:k-2,k-1:k) := U(1:k-2,k:k-1:k) = * = W(1:k-2,kw-1:kw) * ( D(k-1:k,k-1:k)**(-1) ) * IF( K.GT.2 ) THEN * * Compose the columns of the inverse of 2-by-2 pivot * block D in the following way to reduce the number * of FLOPS when we myltiply panel ( W(kw-1) W(kw) ) by * this inverse * * D**(-1) = ( d11 d21 )**(-1) = * ( d21 d22 ) * * = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) = * ( (-d21 ) ( d11 ) ) * * = 1/d21 * 1/((d11/d21)*(d22/d21)-1) * * * * ( ( d22/d21 ) ( -1 ) ) = * ( ( -1 ) ( d11/d21 ) ) * * = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) = * ( ( -1 ) ( D22 ) ) * * = 1/d21 * T * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * * = D21 * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * D21 = W( K-1, KW ) D11 = W( K, KW ) / D21 D22 = W( K-1, KW-1 ) / D21 T = ONE / ( D11*D22-ONE ) D21 = T / D21 * * Update elements in columns A(k-1) and A(k) as * dot products of rows of ( W(kw-1) W(kw) ) and columns * of D**(-1) * DO 20 J = 1, K - 2 A( J, K-1 ) = D21*( D11*W( J, KW-1 )-W( J, KW ) ) A( J, K ) = D21*( D22*W( J, KW )-W( J, KW-1 ) ) 20 CONTINUE END IF * * Copy D(k) to A * A( K-1, K-1 ) = W( K-1, KW-1 ) A( K-1, K ) = W( K-1, KW ) A( K, K ) = W( K, KW ) * END IF * END IF * * Store details of the interchanges in IPIV * IF( KSTEP.EQ.1 ) THEN IPIV( K ) = KP ELSE IPIV( K ) = -KP IPIV( K-1 ) = -KP END IF * * Decrease K and return to the start of the main loop * K = K - KSTEP GO TO 10 * 30 CONTINUE * * Update the upper triangle of A11 (= A(1:k,1:k)) as * * A11 := A11 - U12*D*U12**T = A11 - U12*W**T * * computing blocks of NB columns at a time * DO 50 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB JB = MIN( NB, K-J+1 ) * * Update the upper triangle of the diagonal block * DO 40 JJ = J, J + JB - 1 CALL DGEMV( 'No transpose', JJ-J+1, N-K, -ONE, $ A( J, K+1 ), LDA, W( JJ, KW+1 ), LDW, ONE, $ A( J, JJ ), 1 ) 40 CONTINUE * * Update the rectangular superdiagonal block * CALL DGEMM( 'No transpose', 'Transpose', J-1, JB, N-K, -ONE, $ A( 1, K+1 ), LDA, W( J, KW+1 ), LDW, ONE, $ A( 1, J ), LDA ) 50 CONTINUE * * Put U12 in standard form by partially undoing the interchanges * in columns k+1:n looping backwards from k+1 to n * J = K + 1 60 CONTINUE * * Undo the interchanges (if any) of rows JJ and JP at each * step J * * (Here, J is a diagonal index) JJ = J JP = IPIV( J ) IF( JP.LT.0 ) THEN JP = -JP * (Here, J is a diagonal index) J = J + 1 END IF * (NOTE: Here, J is used to determine row length. Length N-J+1 * of the rows to swap back doesn't include diagonal element) J = J + 1 IF( JP.NE.JJ .AND. J.LE.N ) $ CALL DSWAP( N-J+1, A( JP, J ), LDA, A( JJ, J ), LDA ) IF( J.LT.N ) $ GO TO 60 * * Set KB to the number of columns factorized * KB = N - K * ELSE * * Factorize the leading columns of A using the lower triangle * of A and working forwards, and compute the matrix W = L21*D * for use in updating A22 * * K is the main loop index, increasing from 1 in steps of 1 or 2 * K = 1 70 CONTINUE * * Exit from loop * IF( ( K.GE.NB .AND. NB.LT.N ) .OR. K.GT.N ) $ GO TO 90 * * Copy column K of A to column K of W and update it * CALL DCOPY( N-K+1, A( K, K ), 1, W( K, K ), 1 ) CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ), LDA, $ W( K, 1 ), LDW, ONE, W( K, K ), 1 ) * KSTEP = 1 * * Determine rows and columns to be interchanged and whether * a 1-by-1 or 2-by-2 pivot block will be used * ABSAKK = ABS( W( K, K ) ) * * IMAX is the row-index of the largest off-diagonal element in * column K, and COLMAX is its absolute value. * Determine both COLMAX and IMAX. * IF( K.LT.N ) THEN IMAX = K + IDAMAX( N-K, W( K+1, K ), 1 ) COLMAX = ABS( W( IMAX, K ) ) ELSE COLMAX = ZERO END IF * IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN * * Column K is zero or underflow: set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K KP = K ELSE IF( ABSAKK.GE.ALPHA*COLMAX ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE * * Copy column IMAX to column K+1 of W and update it * CALL DCOPY( IMAX-K, A( IMAX, K ), LDA, W( K, K+1 ), 1 ) CALL DCOPY( N-IMAX+1, A( IMAX, IMAX ), 1, W( IMAX, K+1 ), $ 1 ) CALL DGEMV( 'No transpose', N-K+1, K-1, -ONE, A( K, 1 ), $ LDA, W( IMAX, 1 ), LDW, ONE, W( K, K+1 ), 1 ) * * JMAX is the column-index of the largest off-diagonal * element in row IMAX, and ROWMAX is its absolute value * JMAX = K - 1 + IDAMAX( IMAX-K, W( K, K+1 ), 1 ) ROWMAX = ABS( W( JMAX, K+1 ) ) IF( IMAX.LT.N ) THEN JMAX = IMAX + IDAMAX( N-IMAX, W( IMAX+1, K+1 ), 1 ) ROWMAX = MAX( ROWMAX, ABS( W( JMAX, K+1 ) ) ) END IF * IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE IF( ABS( W( IMAX, K+1 ) ).GE.ALPHA*ROWMAX ) THEN * * interchange rows and columns K and IMAX, use 1-by-1 * pivot block * KP = IMAX * * copy column K+1 of W to column K of W * CALL DCOPY( N-K+1, W( K, K+1 ), 1, W( K, K ), 1 ) ELSE * * interchange rows and columns K+1 and IMAX, use 2-by-2 * pivot block * KP = IMAX KSTEP = 2 END IF END IF * * ============================================================ * * KK is the column of A where pivoting step stopped * KK = K + KSTEP - 1 * * Interchange rows and columns KP and KK. * Updated column KP is already stored in column KK of W. * IF( KP.NE.KK ) THEN * * Copy non-updated column KK to column KP of submatrix A * at step K. No need to copy element into column K * (or K and K+1 for 2-by-2 pivot) of A, since these columns * will be later overwritten. * A( KP, KP ) = A( KK, KK ) CALL DCOPY( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ), $ LDA ) IF( KP.LT.N ) $ CALL DCOPY( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 ) * * Interchange rows KK and KP in first K-1 columns of A * (columns K (or K and K+1 for 2-by-2 pivot) of A will be * later overwritten). Interchange rows KK and KP * in first KK columns of W. * IF( K.GT.1 ) $ CALL DSWAP( K-1, A( KK, 1 ), LDA, A( KP, 1 ), LDA ) CALL DSWAP( KK, W( KK, 1 ), LDW, W( KP, 1 ), LDW ) END IF * IF( KSTEP.EQ.1 ) THEN * * 1-by-1 pivot block D(k): column k of W now holds * * W(k) = L(k)*D(k), * * where L(k) is the k-th column of L * * Store subdiag. elements of column L(k) * and 1-by-1 block D(k) in column k of A. * (NOTE: Diagonal element L(k,k) is a UNIT element * and not stored) * A(k,k) := D(k,k) = W(k,k) * A(k+1:N,k) := L(k+1:N,k) = W(k+1:N,k)/D(k,k) * CALL DCOPY( N-K+1, W( K, K ), 1, A( K, K ), 1 ) IF( K.LT.N ) THEN R1 = ONE / A( K, K ) CALL DSCAL( N-K, R1, A( K+1, K ), 1 ) END IF * ELSE * * 2-by-2 pivot block D(k): columns k and k+1 of W now hold * * ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k) * * where L(k) and L(k+1) are the k-th and (k+1)-th columns * of L * * Store L(k+2:N,k) and L(k+2:N,k+1) and 2-by-2 * block D(k:k+1,k:k+1) in columns k and k+1 of A. * (NOTE: 2-by-2 diagonal block L(k:k+1,k:k+1) is a UNIT * block and not stored) * A(k:k+1,k:k+1) := D(k:k+1,k:k+1) = W(k:k+1,k:k+1) * A(k+2:N,k:k+1) := L(k+2:N,k:k+1) = * = W(k+2:N,k:k+1) * ( D(k:k+1,k:k+1)**(-1) ) * IF( K.LT.N-1 ) THEN * * Compose the columns of the inverse of 2-by-2 pivot * block D in the following way to reduce the number * of FLOPS when we myltiply panel ( W(k) W(k+1) ) by * this inverse * * D**(-1) = ( d11 d21 )**(-1) = * ( d21 d22 ) * * = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) = * ( (-d21 ) ( d11 ) ) * * = 1/d21 * 1/((d11/d21)*(d22/d21)-1) * * * * ( ( d22/d21 ) ( -1 ) ) = * ( ( -1 ) ( d11/d21 ) ) * * = 1/d21 * 1/(D22*D11-1) * ( ( D11 ) ( -1 ) ) = * ( ( -1 ) ( D22 ) ) * * = 1/d21 * T * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * * = D21 * ( ( D11 ) ( -1 ) ) * ( ( -1 ) ( D22 ) ) * D21 = W( K+1, K ) D11 = W( K+1, K+1 ) / D21 D22 = W( K, K ) / D21 T = ONE / ( D11*D22-ONE ) D21 = T / D21 * * Update elements in columns A(k) and A(k+1) as * dot products of rows of ( W(k) W(k+1) ) and columns * of D**(-1) * DO 80 J = K + 2, N A( J, K ) = D21*( D11*W( J, K )-W( J, K+1 ) ) A( J, K+1 ) = D21*( D22*W( J, K+1 )-W( J, K ) ) 80 CONTINUE END IF * * Copy D(k) to A * A( K, K ) = W( K, K ) A( K+1, K ) = W( K+1, K ) A( K+1, K+1 ) = W( K+1, K+1 ) * END IF * END IF * * Store details of the interchanges in IPIV * IF( KSTEP.EQ.1 ) THEN IPIV( K ) = KP ELSE IPIV( K ) = -KP IPIV( K+1 ) = -KP END IF * * Increase K and return to the start of the main loop * K = K + KSTEP GO TO 70 * 90 CONTINUE * * Update the lower triangle of A22 (= A(k:n,k:n)) as * * A22 := A22 - L21*D*L21**T = A22 - L21*W**T * * computing blocks of NB columns at a time * DO 110 J = K, N, NB JB = MIN( NB, N-J+1 ) * * Update the lower triangle of the diagonal block * DO 100 JJ = J, J + JB - 1 CALL DGEMV( 'No transpose', J+JB-JJ, K-1, -ONE, $ A( JJ, 1 ), LDA, W( JJ, 1 ), LDW, ONE, $ A( JJ, JJ ), 1 ) 100 CONTINUE * * Update the rectangular subdiagonal block * IF( J+JB.LE.N ) $ CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB, $ K-1, -ONE, A( J+JB, 1 ), LDA, W( J, 1 ), LDW, $ ONE, A( J+JB, J ), LDA ) 110 CONTINUE * * Put L21 in standard form by partially undoing the interchanges * of rows in columns 1:k-1 looping backwards from k-1 to 1 * J = K - 1 120 CONTINUE * * Undo the interchanges (if any) of rows JJ and JP at each * step J * * (Here, J is a diagonal index) JJ = J JP = IPIV( J ) IF( JP.LT.0 ) THEN JP = -JP * (Here, J is a diagonal index) J = J - 1 END IF * (NOTE: Here, J is used to determine row length. Length J * of the rows to swap back doesn't include diagonal element) J = J - 1 IF( JP.NE.JJ .AND. J.GE.1 ) $ CALL DSWAP( J, A( JP, 1 ), LDA, A( JJ, 1 ), LDA ) IF( J.GT.1 ) $ GO TO 120 * * Set KB to the number of columns factorized * KB = K - 1 * END IF RETURN * * End of DLASYF * END *> \brief \b DSYTF2 computes the factorization of a real symmetric indefinite matrix, using the diagonal pivoting method (unblocked algorithm). * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download DSYTF2 + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE DSYTF2( UPLO, N, A, LDA, IPIV, INFO ) * * .. Scalar Arguments .. * CHARACTER UPLO * INTEGER INFO, LDA, N * .. * .. Array Arguments .. * INTEGER IPIV( * ) * DOUBLE PRECISION A( LDA, * ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DSYTF2 computes the factorization of a real symmetric matrix A using *> the Bunch-Kaufman diagonal pivoting method: *> *> A = U*D*U**T or A = L*D*L**T *> *> where U (or L) is a product of permutation and unit upper (lower) *> triangular matrices, U**T is the transpose of U, and D is symmetric and *> block diagonal with 1-by-1 and 2-by-2 diagonal blocks. *> *> This is the unblocked version of the algorithm, calling Level 2 BLAS. *> \endverbatim * * Arguments: * ========== * *> \param[in] UPLO *> \verbatim *> UPLO is CHARACTER*1 *> Specifies whether the upper or lower triangular part of the *> symmetric matrix A is stored: *> = 'U': Upper triangular *> = 'L': Lower triangular *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> The order of the matrix A. N >= 0. *> \endverbatim *> *> \param[in,out] A *> \verbatim *> A is DOUBLE PRECISION array, dimension (LDA,N) *> On entry, the symmetric matrix A. If UPLO = 'U', the leading *> n-by-n upper triangular part of A contains the upper *> triangular part of the matrix A, and the strictly lower *> triangular part of A is not referenced. If UPLO = 'L', the *> leading n-by-n lower triangular part of A contains the lower *> triangular part of the matrix A, and the strictly upper *> triangular part of A is not referenced. *> *> On exit, the block diagonal matrix D and the multipliers used *> to obtain the factor U or L (see below for further details). *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is INTEGER *> The leading dimension of the array A. LDA >= max(1,N). *> \endverbatim *> *> \param[out] IPIV *> \verbatim *> IPIV is INTEGER array, dimension (N) *> Details of the interchanges and the block structure of D. *> *> If UPLO = 'U': *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were *> interchanged and D(k,k) is a 1-by-1 diagonal block. *> *> If IPIV(k) = IPIV(k-1) < 0, then rows and columns *> k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k) *> is a 2-by-2 diagonal block. *> *> If UPLO = 'L': *> If IPIV(k) > 0, then rows and columns k and IPIV(k) were *> interchanged and D(k,k) is a 1-by-1 diagonal block. *> *> If IPIV(k) = IPIV(k+1) < 0, then rows and columns *> k+1 and -IPIV(k) were interchanged and D(k:k+1,k:k+1) *> is a 2-by-2 diagonal block. *> \endverbatim *> *> \param[out] INFO *> \verbatim *> INFO is INTEGER *> = 0: successful exit *> < 0: if INFO = -k, the k-th argument had an illegal value *> > 0: if INFO = k, D(k,k) is exactly zero. The factorization *> has been completed, but the block diagonal matrix D is *> exactly singular, and division by zero will occur if it *> is used to solve a system of equations. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2013 * *> \ingroup doubleSYcomputational * *> \par Further Details: * ===================== *> *> \verbatim *> *> If UPLO = 'U', then A = U*D*U**T, where *> U = P(n)*U(n)* ... *P(k)U(k)* ..., *> i.e., U is a product of terms P(k)*U(k), where k decreases from n to *> 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as *> defined by IPIV(k), and U(k) is a unit upper triangular matrix, such *> that if the diagonal block D(k) is of order s (s = 1 or 2), then *> *> ( I v 0 ) k-s *> U(k) = ( 0 I 0 ) s *> ( 0 0 I ) n-k *> k-s s n-k *> *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k). *> If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k), *> and A(k,k), and v overwrites A(1:k-2,k-1:k). *> *> If UPLO = 'L', then A = L*D*L**T, where *> L = P(1)*L(1)* ... *P(k)*L(k)* ..., *> i.e., L is a product of terms P(k)*L(k), where k increases from 1 to *> n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1 *> and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as *> defined by IPIV(k), and L(k) is a unit lower triangular matrix, such *> that if the diagonal block D(k) is of order s (s = 1 or 2), then *> *> ( I 0 0 ) k-1 *> L(k) = ( 0 I 0 ) s *> ( 0 v I ) n-k-s+1 *> k-1 s n-k-s+1 *> *> If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k). *> If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k), *> and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1). *> \endverbatim * *> \par Contributors: * ================== *> *> \verbatim *> *> 09-29-06 - patch from *> Bobby Cheng, MathWorks *> *> Replace l.204 and l.372 *> IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN *> by *> IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN *> *> 01-01-96 - Based on modifications by *> J. Lewis, Boeing Computer Services Company *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA *> 1-96 - Based on modifications by J. Lewis, Boeing Computer Services *> Company *> \endverbatim * * ===================================================================== SUBROUTINE DSYTF2( UPLO, N, A, LDA, IPIV, INFO ) * * -- LAPACK computational routine (version 3.5.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2013 * * .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, LDA, N * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ) * .. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) DOUBLE PRECISION EIGHT, SEVTEN PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 ) * .. * .. Local Scalars .. LOGICAL UPPER INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1, $ ROWMAX, T, WK, WKM1, WKP1 * .. * .. External Functions .. LOGICAL LSAME, DISNAN INTEGER IDAMAX EXTERNAL LSAME, IDAMAX, DISNAN * .. * .. External Subroutines .. EXTERNAL DSCAL, DSWAP, DSYR, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DSYTF2', -INFO ) RETURN END IF * * Initialize ALPHA for use in choosing pivot block size. * ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT * IF( UPPER ) THEN * * Factorize A as U*D*U**T using the upper triangle of A * * K is the main loop index, decreasing from N to 1 in steps of * 1 or 2 * K = N 10 CONTINUE * * If K < 1, exit from loop * IF( K.LT.1 ) $ GO TO 70 KSTEP = 1 * * Determine rows and columns to be interchanged and whether * a 1-by-1 or 2-by-2 pivot block will be used * ABSAKK = ABS( A( K, K ) ) * * IMAX is the row-index of the largest off-diagonal element in * column K, and COLMAX is its absolute value. * Determine both COLMAX and IMAX. * IF( K.GT.1 ) THEN IMAX = IDAMAX( K-1, A( 1, K ), 1 ) COLMAX = ABS( A( IMAX, K ) ) ELSE COLMAX = ZERO END IF * IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN * * Column K is zero or underflow, or contains a NaN: * set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K KP = K ELSE IF( ABSAKK.GE.ALPHA*COLMAX ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE * * JMAX is the column-index of the largest off-diagonal * element in row IMAX, and ROWMAX is its absolute value * JMAX = IMAX + IDAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA ) ROWMAX = ABS( A( IMAX, JMAX ) ) IF( IMAX.GT.1 ) THEN JMAX = IDAMAX( IMAX-1, A( 1, IMAX ), 1 ) ROWMAX = MAX( ROWMAX, ABS( A( JMAX, IMAX ) ) ) END IF * IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE IF( ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN * * interchange rows and columns K and IMAX, use 1-by-1 * pivot block * KP = IMAX ELSE * * interchange rows and columns K-1 and IMAX, use 2-by-2 * pivot block * KP = IMAX KSTEP = 2 END IF END IF * KK = K - KSTEP + 1 IF( KP.NE.KK ) THEN * * Interchange rows and columns KK and KP in the leading * submatrix A(1:k,1:k) * CALL DSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 ) CALL DSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ), $ LDA ) T = A( KK, KK ) A( KK, KK ) = A( KP, KP ) A( KP, KP ) = T IF( KSTEP.EQ.2 ) THEN T = A( K-1, K ) A( K-1, K ) = A( KP, K ) A( KP, K ) = T END IF END IF * * Update the leading submatrix * IF( KSTEP.EQ.1 ) THEN * * 1-by-1 pivot block D(k): column k now holds * * W(k) = U(k)*D(k) * * where U(k) is the k-th column of U * * Perform a rank-1 update of A(1:k-1,1:k-1) as * * A := A - U(k)*D(k)*U(k)**T = A - W(k)*1/D(k)*W(k)**T * R1 = ONE / A( K, K ) CALL DSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA ) * * Store U(k) in column k * CALL DSCAL( K-1, R1, A( 1, K ), 1 ) ELSE * * 2-by-2 pivot block D(k): columns k and k-1 now hold * * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k) * * where U(k) and U(k-1) are the k-th and (k-1)-th columns * of U * * Perform a rank-2 update of A(1:k-2,1:k-2) as * * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )**T * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )**T * IF( K.GT.2 ) THEN * D12 = A( K-1, K ) D22 = A( K-1, K-1 ) / D12 D11 = A( K, K ) / D12 T = ONE / ( D11*D22-ONE ) D12 = T / D12 * DO 30 J = K - 2, 1, -1 WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) ) WK = D12*( D22*A( J, K )-A( J, K-1 ) ) DO 20 I = J, 1, -1 A( I, J ) = A( I, J ) - A( I, K )*WK - $ A( I, K-1 )*WKM1 20 CONTINUE A( J, K ) = WK A( J, K-1 ) = WKM1 30 CONTINUE * END IF * END IF END IF * * Store details of the interchanges in IPIV * IF( KSTEP.EQ.1 ) THEN IPIV( K ) = KP ELSE IPIV( K ) = -KP IPIV( K-1 ) = -KP END IF * * Decrease K and return to the start of the main loop * K = K - KSTEP GO TO 10 * ELSE * * Factorize A as L*D*L**T using the lower triangle of A * * K is the main loop index, increasing from 1 to N in steps of * 1 or 2 * K = 1 40 CONTINUE * * If K > N, exit from loop * IF( K.GT.N ) $ GO TO 70 KSTEP = 1 * * Determine rows and columns to be interchanged and whether * a 1-by-1 or 2-by-2 pivot block will be used * ABSAKK = ABS( A( K, K ) ) * * IMAX is the row-index of the largest off-diagonal element in * column K, and COLMAX is its absolute value. * Determine both COLMAX and IMAX. * IF( K.LT.N ) THEN IMAX = K + IDAMAX( N-K, A( K+1, K ), 1 ) COLMAX = ABS( A( IMAX, K ) ) ELSE COLMAX = ZERO END IF * IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN * * Column K is zero or underflow, or contains a NaN: * set INFO and continue * IF( INFO.EQ.0 ) $ INFO = K KP = K ELSE IF( ABSAKK.GE.ALPHA*COLMAX ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE * * JMAX is the column-index of the largest off-diagonal * element in row IMAX, and ROWMAX is its absolute value * JMAX = K - 1 + IDAMAX( IMAX-K, A( IMAX, K ), LDA ) ROWMAX = ABS( A( IMAX, JMAX ) ) IF( IMAX.LT.N ) THEN JMAX = IMAX + IDAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 ) ROWMAX = MAX( ROWMAX, ABS( A( JMAX, IMAX ) ) ) END IF * IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN * * no interchange, use 1-by-1 pivot block * KP = K ELSE IF( ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN * * interchange rows and columns K and IMAX, use 1-by-1 * pivot block * KP = IMAX ELSE * * interchange rows and columns K+1 and IMAX, use 2-by-2 * pivot block * KP = IMAX KSTEP = 2 END IF END IF * KK = K + KSTEP - 1 IF( KP.NE.KK ) THEN * * Interchange rows and columns KK and KP in the trailing * submatrix A(k:n,k:n) * IF( KP.LT.N ) $ CALL DSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 ) CALL DSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ), $ LDA ) T = A( KK, KK ) A( KK, KK ) = A( KP, KP ) A( KP, KP ) = T IF( KSTEP.EQ.2 ) THEN T = A( K+1, K ) A( K+1, K ) = A( KP, K ) A( KP, K ) = T END IF END IF * * Update the trailing submatrix * IF( KSTEP.EQ.1 ) THEN * * 1-by-1 pivot block D(k): column k now holds * * W(k) = L(k)*D(k) * * where L(k) is the k-th column of L * IF( K.LT.N ) THEN * * Perform a rank-1 update of A(k+1:n,k+1:n) as * * A := A - L(k)*D(k)*L(k)**T = A - W(k)*(1/D(k))*W(k)**T * D11 = ONE / A( K, K ) CALL DSYR( UPLO, N-K, -D11, A( K+1, K ), 1, $ A( K+1, K+1 ), LDA ) * * Store L(k) in column K * CALL DSCAL( N-K, D11, A( K+1, K ), 1 ) END IF ELSE * * 2-by-2 pivot block D(k) * IF( K.LT.N-1 ) THEN * * Perform a rank-2 update of A(k+2:n,k+2:n) as * * A := A - ( (A(k) A(k+1))*D(k)**(-1) ) * (A(k) A(k+1))**T * * where L(k) and L(k+1) are the k-th and (k+1)-th * columns of L * D21 = A( K+1, K ) D11 = A( K+1, K+1 ) / D21 D22 = A( K, K ) / D21 T = ONE / ( D11*D22-ONE ) D21 = T / D21 * DO 60 J = K + 2, N * WK = D21*( D11*A( J, K )-A( J, K+1 ) ) WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) ) * DO 50 I = J, N A( I, J ) = A( I, J ) - A( I, K )*WK - $ A( I, K+1 )*WKP1 50 CONTINUE * A( J, K ) = WK A( J, K+1 ) = WKP1 * 60 CONTINUE END IF END IF END IF * * Store details of the interchanges in IPIV * IF( KSTEP.EQ.1 ) THEN IPIV( K ) = KP ELSE IPIV( K ) = -KP IPIV( K+1 ) = -KP END IF * * Increase K and return to the start of the main loop * K = K + KSTEP GO TO 40 * END IF * 70 CONTINUE * RETURN * * End of DSYTF2 * END *> \brief \b XERBLA * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download XERBLA + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * SUBROUTINE XERBLA( SRNAME, INFO ) * * .. Scalar Arguments .. * CHARACTER*(*) SRNAME * INTEGER INFO * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> 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. *> \endverbatim * * Arguments: * ========== * *> \param[in] SRNAME *> \verbatim *> SRNAME is CHARACTER*(*) *> The name of the routine which called XERBLA. *> \endverbatim *> *> \param[in] INFO *> \verbatim *> INFO is INTEGER *> The position of the invalid parameter in the parameter list *> of the calling routine. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2011 * *> \ingroup auxOTHERauxiliary * * ===================================================================== SUBROUTINE XERBLA( SRNAME, INFO ) * * -- LAPACK auxiliary routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Scalar Arguments .. CHARACTER*(*) SRNAME INTEGER INFO * .. * * ===================================================================== * * .. Intrinsic Functions .. INTRINSIC LEN_TRIM * .. * .. Executable Statements .. * WRITE( *, FMT = 9999 )SRNAME( 1:LEN_TRIM( SRNAME ) ), INFO * STOP * 9999 FORMAT( ' ** On entry to ', A, ' parameter number ', I2, ' had ', $ 'an illegal value' ) * * End of XERBLA * END *> \brief \b ILAENV * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download ILAENV + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, N4 ) * * .. Scalar Arguments .. * CHARACTER*( * ) NAME, OPTS * INTEGER ISPEC, N1, N2, N3, N4 * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> ILAENV is called from the LAPACK routines to choose problem-dependent *> parameters for the local environment. See ISPEC for a description of *> the parameters. *> *> ILAENV returns an INTEGER *> if ILAENV >= 0: ILAENV returns the value of the parameter specified by ISPEC *> if ILAENV < 0: if ILAENV = -k, the k-th argument had an illegal value. *> *> 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. *> \endverbatim * * Arguments: * ========== * *> \param[in] ISPEC *> \verbatim *> ISPEC is 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 (DEPRECATED) *> = 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 method *> for nonsymmetric eigenvalue problems (DEPRECATED) *> = 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 *> 12 <= ISPEC <= 16: *> xHSEQR or related subroutines, *> see IPARMQ for detailed explanation *> \endverbatim *> *> \param[in] NAME *> \verbatim *> NAME is CHARACTER*(*) *> The name of the calling subroutine, in either upper case or *> lower case. *> \endverbatim *> *> \param[in] OPTS *> \verbatim *> OPTS is 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'. *> \endverbatim *> *> \param[in] N1 *> \verbatim *> N1 is INTEGER *> \endverbatim *> *> \param[in] N2 *> \verbatim *> N2 is INTEGER *> \endverbatim *> *> \param[in] N3 *> \verbatim *> N3 is INTEGER *> \endverbatim *> *> \param[in] N4 *> \verbatim *> N4 is INTEGER *> Problem dimensions for the subroutine NAME; these may not all *> be required. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date June 2016 * *> \ingroup auxOTHERauxiliary * *> \par Further Details: * ===================== *> *> \verbatim *> *> 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 ) *> \endverbatim *> * ===================================================================== INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, N4 ) * * -- LAPACK auxiliary routine (version 3.6.1) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * June 2016 * * .. Scalar Arguments .. CHARACTER*( * ) NAME, OPTS INTEGER ISPEC, N1, N2, N3, N4 * .. * * ===================================================================== * * .. Local Scalars .. INTEGER I, IC, IZ, NB, NBMIN, NX LOGICAL CNAME, SNAME CHARACTER C1*1, C2*2, C4*2, C3*3, SUBNAM*6 * .. * .. Intrinsic Functions .. INTRINSIC CHAR, ICHAR, INT, MIN, REAL * .. * .. External Functions .. INTEGER IEEECK, IPARMQ EXTERNAL IEEECK, IPARMQ * .. * .. Executable Statements .. * GO TO ( 10, 10, 10, 80, 90, 100, 110, 120, $ 130, 140, 150, 160, 160, 160, 160, 160 )ISPEC * * Invalid value for ISPEC * ILAENV = -1 RETURN * 10 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 20 I = 2, 6 IC = ICHAR( SUBNAM( I: I ) ) IF( IC.GE.97 .AND. IC.LE.122 ) $ SUBNAM( I: I ) = CHAR( IC-32 ) 20 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 30 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 ) 30 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 40 I = 2, 6 IC = ICHAR( SUBNAM( I: I ) ) IF( IC.GE.225 .AND. IC.LE.250 ) $ SUBNAM( I: I ) = CHAR( IC-32 ) 40 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 ( 50, 60, 70 )ISPEC * 50 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 double precision. * 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 ELSE IF ( C3.EQ.'EVC' ) 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 ELSE IF( C2.EQ.'GG' ) THEN NB = 32 IF( C3.EQ.'HD3' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF END IF END IF ILAENV = NB RETURN * 60 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 ELSE IF( C2.EQ.'GG' ) THEN NBMIN = 2 IF( C3.EQ.'HD3' ) THEN NBMIN = 2 END IF END IF ILAENV = NBMIN RETURN * 70 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 ELSE IF( C2.EQ.'GG' ) THEN NX = 128 IF( C3.EQ.'HD3' ) THEN NX = 128 END IF END IF ILAENV = NX RETURN * 80 CONTINUE * * ISPEC = 4: number of shifts (used by xHSEQR) * ILAENV = 6 RETURN * 90 CONTINUE * * ISPEC = 5: minimum column dimension (not used) * ILAENV = 2 RETURN * 100 CONTINUE * * ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) * ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 ) RETURN * 110 CONTINUE * * ISPEC = 7: number of processors (not used) * ILAENV = 1 RETURN * 120 CONTINUE * * ISPEC = 8: crossover point for multishift (used by xHSEQR) * ILAENV = 50 RETURN * 130 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 * 140 CONTINUE * * ISPEC = 10: ieee NaN arithmetic can be trusted not to trap * * ILAENV = 0 ILAENV = 1 IF( ILAENV.EQ.1 ) THEN ILAENV = IEEECK( 1, 0.0, 1.0 ) END IF RETURN * 150 CONTINUE * * ISPEC = 11: infinity arithmetic can be trusted not to trap * * ILAENV = 0 ILAENV = 1 IF( ILAENV.EQ.1 ) THEN ILAENV = IEEECK( 0, 0.0, 1.0 ) END IF RETURN * 160 CONTINUE * * 12 <= ISPEC <= 16: xHSEQR or related subroutines. * ILAENV = IPARMQ( ISPEC, NAME, OPTS, N1, N2, N3, N4 ) RETURN * * End of ILAENV * END *> \brief \b DSWAP * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * * Definition: * =========== * * SUBROUTINE DSWAP(N,DX,INCX,DY,INCY) * * .. Scalar Arguments .. * INTEGER INCX,INCY,N * .. * .. Array Arguments .. * DOUBLE PRECISION DX(*),DY(*) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> interchanges two vectors. *> uses unrolled loops for increments equal one. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2011 * *> \ingroup double_blas_level1 * *> \par Further Details: * ===================== *> *> \verbatim *> *> jack dongarra, linpack, 3/11/78. *> modified 12/3/93, array(1) declarations changed to array(*) *> \endverbatim *> * ===================================================================== SUBROUTINE DSWAP(N,DX,INCX,DY,INCY) * * -- Reference BLAS level1 routine (version 3.4.0) -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Scalar Arguments .. INTEGER INCX,INCY,N * .. * .. Array Arguments .. DOUBLE PRECISION DX(*),DY(*) * .. * * ===================================================================== * * .. Local Scalars .. DOUBLE PRECISION DTEMP INTEGER I,IX,IY,M,MP1 * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. IF (N.LE.0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN * * code for both increments equal to 1 * * * clean-up loop * M = MOD(N,3) IF (M.NE.0) THEN DO I = 1,M DTEMP = DX(I) DX(I) = DY(I) DY(I) = DTEMP END DO IF (N.LT.3) RETURN END IF MP1 = M + 1 DO 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 END DO ELSE * * code for unequal increments or equal increments not equal * to 1 * IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO I = 1,N DTEMP = DX(IX) DX(IX) = DY(IY) DY(IY) = DTEMP IX = IX + INCX IY = IY + INCY END DO END IF RETURN END *> \brief \b DCOPY * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * * Definition: * =========== * * SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) * * .. Scalar Arguments .. * INTEGER INCX,INCY,N * .. * .. Array Arguments .. * DOUBLE PRECISION DX(*),DY(*) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DCOPY copies a vector, x, to a vector, y. *> uses unrolled loops for increments equal to one. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2011 * *> \ingroup double_blas_level1 * *> \par Further Details: * ===================== *> *> \verbatim *> *> jack dongarra, linpack, 3/11/78. *> modified 12/3/93, array(1) declarations changed to array(*) *> \endverbatim *> * ===================================================================== SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) * * -- Reference BLAS level1 routine (version 3.4.0) -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Scalar Arguments .. INTEGER INCX,INCY,N * .. * .. Array Arguments .. DOUBLE PRECISION DX(*),DY(*) * .. * * ===================================================================== * * .. Local Scalars .. INTEGER I,IX,IY,M,MP1 * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. IF (N.LE.0) RETURN IF (INCX.EQ.1 .AND. INCY.EQ.1) THEN * * code for both increments equal to 1 * * * clean-up loop * M = MOD(N,7) IF (M.NE.0) THEN DO I = 1,M DY(I) = DX(I) END DO IF (N.LT.7) RETURN END IF MP1 = M + 1 DO I = MP1,N,7 DY(I) = DX(I) DY(I+1) = DX(I+1) DY(I+2) = DX(I+2) DY(I+3) = DX(I+3) DY(I+4) = DX(I+4) DY(I+5) = DX(I+5) DY(I+6) = DX(I+6) END DO ELSE * * code for unequal increments or equal increments * not equal to 1 * IX = 1 IY = 1 IF (INCX.LT.0) IX = (-N+1)*INCX + 1 IF (INCY.LT.0) IY = (-N+1)*INCY + 1 DO I = 1,N DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY END DO END IF RETURN END *> \brief \b DSCAL * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * * Definition: * =========== * * SUBROUTINE DSCAL(N,DA,DX,INCX) * * .. Scalar Arguments .. * DOUBLE PRECISION DA * INTEGER INCX,N * .. * .. Array Arguments .. * DOUBLE PRECISION DX(*) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DSCAL scales a vector by a constant. *> uses unrolled loops for increment equal to one. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2011 * *> \ingroup double_blas_level1 * *> \par Further Details: * ===================== *> *> \verbatim *> *> jack dongarra, linpack, 3/11/78. *> modified 3/93 to return if incx .le. 0. *> modified 12/3/93, array(1) declarations changed to array(*) *> \endverbatim *> * ===================================================================== SUBROUTINE DSCAL(N,DA,DX,INCX) * * -- Reference BLAS level1 routine (version 3.4.0) -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Scalar Arguments .. DOUBLE PRECISION DA INTEGER INCX,N * .. * .. Array Arguments .. DOUBLE PRECISION DX(*) * .. * * ===================================================================== * * .. Local Scalars .. INTEGER I,M,MP1,NINCX * .. * .. Intrinsic Functions .. INTRINSIC MOD * .. IF (N.LE.0 .OR. INCX.LE.0) RETURN IF (INCX.EQ.1) THEN * * code for increment equal to 1 * * * clean-up loop * M = MOD(N,5) IF (M.NE.0) THEN DO I = 1,M DX(I) = DA*DX(I) END DO IF (N.LT.5) RETURN END IF MP1 = M + 1 DO 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) END DO ELSE * * code for increment not equal to 1 * NINCX = N*INCX DO I = 1,NINCX,INCX DX(I) = DA*DX(I) END DO END IF RETURN END *> \brief \b IDAMAX * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * * Definition: * =========== * * INTEGER FUNCTION IDAMAX(N,DX,INCX) * * .. Scalar Arguments .. * INTEGER INCX,N * .. * .. Array Arguments .. * DOUBLE PRECISION DX(*) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> IDAMAX finds the index of the first element having maximum absolute value. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2015 * *> \ingroup aux_blas * *> \par Further Details: * ===================== *> *> \verbatim *> *> jack dongarra, linpack, 3/11/78. *> modified 3/93 to return if incx .le. 0. *> modified 12/3/93, array(1) declarations changed to array(*) *> \endverbatim *> * ===================================================================== INTEGER FUNCTION IDAMAX(N,DX,INCX) * * -- Reference BLAS level1 routine (version 3.6.0) -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2015 * * .. Scalar Arguments .. INTEGER INCX,N * .. * .. Array Arguments .. DOUBLE PRECISION DX(*) * .. * * ===================================================================== * * .. Local Scalars .. DOUBLE PRECISION DMAX INTEGER I,IX * .. * .. Intrinsic Functions .. INTRINSIC DABS * .. IDAMAX = 0 IF (N.LT.1 .OR. INCX.LE.0) RETURN IDAMAX = 1 IF (N.EQ.1) RETURN IF (INCX.EQ.1) THEN * * code for increment equal to 1 * DMAX = DABS(DX(1)) DO I = 2,N IF (DABS(DX(I)).GT.DMAX) THEN IDAMAX = I DMAX = DABS(DX(I)) END IF END DO ELSE * * code for increment not equal to 1 * IX = 1 DMAX = DABS(DX(1)) IX = IX + INCX DO I = 2,N IF (DABS(DX(IX)).GT.DMAX) THEN IDAMAX = I DMAX = DABS(DX(IX)) END IF IX = IX + INCX END DO END IF RETURN END *> \brief \b IPARMQ * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download IPARMQ + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * INTEGER FUNCTION IPARMQ( ISPEC, NAME, OPTS, N, ILO, IHI, LWORK ) * * .. Scalar Arguments .. * INTEGER IHI, ILO, ISPEC, LWORK, N * CHARACTER NAME*( * ), OPTS*( * ) * * *> \par Purpose: * ============= *> *> \verbatim *> *> This program sets problem and machine dependent parameters *> useful for xHSEQR and related subroutines for eigenvalue *> problems. It is called whenever *> IPARMQ is called with 12 <= ISPEC <= 16 *> \endverbatim * * Arguments: * ========== * *> \param[in] ISPEC *> \verbatim *> ISPEC is integer scalar *> ISPEC specifies which tunable parameter IPARMQ should *> return. *> *> ISPEC=12: (INMIN) Matrices of order nmin or less *> are sent directly to xLAHQR, the implicit *> double shift QR algorithm. NMIN must be *> at least 11. *> *> ISPEC=13: (INWIN) Size of the deflation window. *> This is best set greater than or equal to *> the number of simultaneous shifts NS. *> Larger matrices benefit from larger deflation *> windows. *> *> ISPEC=14: (INIBL) Determines when to stop nibbling and *> invest in an (expensive) multi-shift QR sweep. *> If the aggressive early deflation subroutine *> finds LD converged eigenvalues from an order *> NW deflation window and LD.GT.(NW*NIBBLE)/100, *> then the next QR sweep is skipped and early *> deflation is applied immediately to the *> remaining active diagonal block. Setting *> IPARMQ(ISPEC=14) = 0 causes TTQRE to skip a *> multi-shift QR sweep whenever early deflation *> finds a converged eigenvalue. Setting *> IPARMQ(ISPEC=14) greater than or equal to 100 *> prevents TTQRE from skipping a multi-shift *> QR sweep. *> *> ISPEC=15: (NSHFTS) The number of simultaneous shifts in *> a multi-shift QR iteration. *> *> ISPEC=16: (IACC22) IPARMQ is set to 0, 1 or 2 with the *> following meanings. *> 0: During the multi-shift QR/QZ sweep, *> blocked eigenvalue reordering, blocked *> Hessenberg-triangular reduction, *> reflections and/or rotations are not *> accumulated when updating the *> far-from-diagonal matrix entries. *> 1: During the multi-shift QR/QZ sweep, *> blocked eigenvalue reordering, blocked *> Hessenberg-triangular reduction, *> reflections and/or rotations are *> accumulated, and matrix-matrix *> multiplication is used to update the *> far-from-diagonal matrix entries. *> 2: During the multi-shift QR/QZ sweep, *> blocked eigenvalue reordering, blocked *> Hessenberg-triangular reduction, *> reflections and/or rotations are *> accumulated, and 2-by-2 block structure *> is exploited during matrix-matrix *> multiplies. *> (If xTRMM is slower than xGEMM, then *> IPARMQ(ISPEC=16)=1 may be more efficient than *> IPARMQ(ISPEC=16)=2 despite the greater level of *> arithmetic work implied by the latter choice.) *> \endverbatim *> *> \param[in] NAME *> \verbatim *> NAME is character string *> Name of the calling subroutine *> \endverbatim *> *> \param[in] OPTS *> \verbatim *> OPTS is character string *> This is a concatenation of the string arguments to *> TTQRE. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is integer scalar *> N is the order of the Hessenberg matrix H. *> \endverbatim *> *> \param[in] ILO *> \verbatim *> ILO is INTEGER *> \endverbatim *> *> \param[in] IHI *> \verbatim *> IHI is INTEGER *> It is assumed that H is already upper triangular *> in rows and columns 1:ILO-1 and IHI+1:N. *> \endverbatim *> *> \param[in] LWORK *> \verbatim *> LWORK is integer scalar *> The amount of workspace available. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2015 * *> \ingroup auxOTHERauxiliary * *> \par Further Details: * ===================== *> *> \verbatim *> *> Little is known about how best to choose these parameters. *> It is possible to use different values of the parameters *> for each of CHSEQR, DHSEQR, SHSEQR and ZHSEQR. *> *> It is probably best to choose different parameters for *> different matrices and different parameters at different *> times during the iteration, but this has not been *> implemented --- yet. *> *> *> The best choices of most of the parameters depend *> in an ill-understood way on the relative execution *> rate of xLAQR3 and xLAQR5 and on the nature of each *> particular eigenvalue problem. Experiment may be the *> only practical way to determine which choices are most *> effective. *> *> Following is a list of default values supplied by IPARMQ. *> These defaults may be adjusted in order to attain better *> performance in any particular computational environment. *> *> IPARMQ(ISPEC=12) The xLAHQR vs xLAQR0 crossover point. *> Default: 75. (Must be at least 11.) *> *> IPARMQ(ISPEC=13) Recommended deflation window size. *> This depends on ILO, IHI and NS, the *> number of simultaneous shifts returned *> by IPARMQ(ISPEC=15). The default for *> (IHI-ILO+1).LE.500 is NS. The default *> for (IHI-ILO+1).GT.500 is 3*NS/2. *> *> IPARMQ(ISPEC=14) Nibble crossover point. Default: 14. *> *> IPARMQ(ISPEC=15) Number of simultaneous shifts, NS. *> a multi-shift QR iteration. *> *> If IHI-ILO+1 is ... *> *> greater than ...but less ... the *> or equal to ... than default is *> *> 0 30 NS = 2+ *> 30 60 NS = 4+ *> 60 150 NS = 10 *> 150 590 NS = ** *> 590 3000 NS = 64 *> 3000 6000 NS = 128 *> 6000 infinity NS = 256 *> *> (+) By default matrices of this order are *> passed to the implicit double shift routine *> xLAHQR. See IPARMQ(ISPEC=12) above. These *> values of NS are used only in case of a rare *> xLAHQR failure. *> *> (**) The asterisks (**) indicate an ad-hoc *> function increasing from 10 to 64. *> *> IPARMQ(ISPEC=16) Select structured matrix multiply. *> (See ISPEC=16 above for details.) *> Default: 3. *> \endverbatim *> * ===================================================================== INTEGER FUNCTION IPARMQ( ISPEC, NAME, OPTS, N, ILO, IHI, LWORK ) * * -- LAPACK auxiliary routine (version 3.6.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2015 * * .. Scalar Arguments .. INTEGER IHI, ILO, ISPEC, LWORK, N CHARACTER NAME*( * ), OPTS*( * ) * * ================================================================ * .. Parameters .. INTEGER INMIN, INWIN, INIBL, ISHFTS, IACC22 PARAMETER ( INMIN = 12, INWIN = 13, INIBL = 14, $ ISHFTS = 15, IACC22 = 16 ) INTEGER NMIN, K22MIN, KACMIN, NIBBLE, KNWSWP PARAMETER ( NMIN = 75, K22MIN = 14, KACMIN = 14, $ NIBBLE = 14, KNWSWP = 500 ) REAL TWO PARAMETER ( TWO = 2.0 ) * .. * .. Local Scalars .. INTEGER NH, NS INTEGER I, IC, IZ CHARACTER SUBNAM*6 * .. * .. Intrinsic Functions .. INTRINSIC LOG, MAX, MOD, NINT, REAL * .. * .. Executable Statements .. IF( ( ISPEC.EQ.ISHFTS ) .OR. ( ISPEC.EQ.INWIN ) .OR. $ ( ISPEC.EQ.IACC22 ) ) THEN * * ==== Set the number simultaneous shifts ==== * NH = IHI - ILO + 1 NS = 2 IF( NH.GE.30 ) $ NS = 4 IF( NH.GE.60 ) $ NS = 10 IF( NH.GE.150 ) $ NS = MAX( 10, NH / NINT( LOG( REAL( NH ) ) / LOG( TWO ) ) ) IF( NH.GE.590 ) $ NS = 64 IF( NH.GE.3000 ) $ NS = 128 IF( NH.GE.6000 ) $ NS = 256 NS = MAX( 2, NS-MOD( NS, 2 ) ) END IF * IF( ISPEC.EQ.INMIN ) THEN * * * ===== Matrices of order smaller than NMIN get sent * . to xLAHQR, the classic double shift algorithm. * . This must be at least 11. ==== * IPARMQ = NMIN * ELSE IF( ISPEC.EQ.INIBL ) THEN * * ==== INIBL: skip a multi-shift qr iteration and * . whenever aggressive early deflation finds * . at least (NIBBLE*(window size)/100) deflations. ==== * IPARMQ = NIBBLE * ELSE IF( ISPEC.EQ.ISHFTS ) THEN * * ==== NSHFTS: The number of simultaneous shifts ===== * IPARMQ = NS * ELSE IF( ISPEC.EQ.INWIN ) THEN * * ==== NW: deflation window size. ==== * IF( NH.LE.KNWSWP ) THEN IPARMQ = NS ELSE IPARMQ = 3*NS / 2 END IF * ELSE IF( ISPEC.EQ.IACC22 ) THEN * * ==== IACC22: Whether to accumulate reflections * . before updating the far-from-diagonal elements * . and whether to use 2-by-2 block structure while * . doing it. A small amount of work could be saved * . by making this choice dependent also upon the * . NH=IHI-ILO+1. * * * Convert NAME to upper case if the first character is lower case. * IPARMQ = 0 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 I = 2, 6 IC = ICHAR( SUBNAM( I: I ) ) IF( IC.GE.97 .AND. IC.LE.122 ) $ SUBNAM( I: I ) = CHAR( IC-32 ) END DO 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 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 ) END DO 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 I = 2, 6 IC = ICHAR( SUBNAM( I: I ) ) IF( IC.GE.225 .AND. IC.LE.250 ) $ SUBNAM( I: I ) = CHAR( IC-32 ) END DO END IF END IF * IF( SUBNAM( 2:6 ).EQ.'GGHRD' .OR. $ SUBNAM( 2:6 ).EQ.'GGHD3' ) THEN IPARMQ = 1 IF( NH.GE.K22MIN ) $ IPARMQ = 2 ELSE IF ( SUBNAM( 4:6 ).EQ.'EXC' ) THEN IF( NH.GE.KACMIN ) $ IPARMQ = 1 IF( NH.GE.K22MIN ) $ IPARMQ = 2 ELSE IF ( SUBNAM( 2:6 ).EQ.'HSEQR' .OR. $ SUBNAM( 2:5 ).EQ.'LAQR' ) THEN IF( NS.GE.KACMIN ) $ IPARMQ = 1 IF( NS.GE.K22MIN ) $ IPARMQ = 2 END IF * ELSE * ===== invalid value of ispec ===== IPARMQ = -1 * END IF * * ==== End of IPARMQ ==== * END *> \brief \b IEEECK * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download IEEECK + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * INTEGER FUNCTION IEEECK( ISPEC, ZERO, ONE ) * * .. Scalar Arguments .. * INTEGER ISPEC * REAL ONE, ZERO * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> IEEECK is called from the ILAENV to verify that Infinity and *> possibly NaN arithmetic is safe (i.e. will not trap). *> \endverbatim * * Arguments: * ========== * *> \param[in] ISPEC *> \verbatim *> ISPEC is 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. *> \endverbatim *> *> \param[in] ZERO *> \verbatim *> ZERO is REAL *> Must contain the value 0.0 *> This is passed to prevent the compiler from optimizing *> away this code. *> \endverbatim *> *> \param[in] ONE *> \verbatim *> ONE is 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 *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2011 * *> \ingroup auxOTHERauxiliary * * ===================================================================== INTEGER FUNCTION IEEECK( ISPEC, ZERO, ONE ) * * -- LAPACK auxiliary routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Scalar Arguments .. INTEGER ISPEC REAL ONE, ZERO * .. * * ===================================================================== * * .. Local Scalars .. REAL 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*ZERO * 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 *> \brief \b DGEMV * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * * Definition: * =========== * * SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) * * .. Scalar Arguments .. * DOUBLE PRECISION ALPHA,BETA * INTEGER INCX,INCY,LDA,M,N * CHARACTER TRANS * .. * .. Array Arguments .. * DOUBLE PRECISION A(LDA,*),X(*),Y(*) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DGEMV performs one of the matrix-vector operations *> *> y := alpha*A*x + beta*y, or y := alpha*A**T*x + beta*y, *> *> where alpha and beta are scalars, x and y are vectors and A is an *> m by n matrix. *> \endverbatim * * Arguments: * ========== * *> \param[in] TRANS *> \verbatim *> TRANS is 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**T*x + beta*y. *> *> TRANS = 'C' or 'c' y := alpha*A**T*x + beta*y. *> \endverbatim *> *> \param[in] M *> \verbatim *> M is INTEGER *> On entry, M specifies the number of rows of the matrix A. *> M must be at least zero. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> On entry, N specifies the number of columns of the matrix A. *> N must be at least zero. *> \endverbatim *> *> \param[in] ALPHA *> \verbatim *> ALPHA is DOUBLE PRECISION. *> On entry, ALPHA specifies the scalar alpha. *> \endverbatim *> *> \param[in] A *> \verbatim *> A is 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. *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is 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 ). *> \endverbatim *> *> \param[in] X *> \verbatim *> X is 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. *> \endverbatim *> *> \param[in] INCX *> \verbatim *> INCX is INTEGER *> On entry, INCX specifies the increment for the elements of *> X. INCX must not be zero. *> \endverbatim *> *> \param[in] BETA *> \verbatim *> BETA is DOUBLE PRECISION. *> On entry, BETA specifies the scalar beta. When BETA is *> supplied as zero then Y need not be set on input. *> \endverbatim *> *> \param[in,out] Y *> \verbatim *> Y is 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. *> \endverbatim *> *> \param[in] INCY *> \verbatim *> INCY is INTEGER *> On entry, INCY specifies the increment for the elements of *> Y. INCY must not be zero. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2015 * *> \ingroup double_blas_level2 * *> \par Further Details: * ===================== *> *> \verbatim *> *> Level 2 Blas routine. *> The vector and matrix arguments are not referenced when N = 0, or M = 0 *> *> -- 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. *> \endverbatim *> * ===================================================================== SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) * * -- Reference BLAS level2 routine (version 3.6.0) -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2015 * * .. Scalar Arguments .. DOUBLE PRECISION ALPHA,BETA INTEGER INCX,INCY,LDA,M,N CHARACTER TRANS * .. * .. Array Arguments .. DOUBLE PRECISION A(LDA,*),X(*),Y(*) * .. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE,ZERO PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) * .. * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY,LENX,LENY * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * * Test the input parameters. * INFO = 0 IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND. + .NOT.LSAME(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 (LSAME(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 (LSAME(TRANS,'N')) THEN * * Form y := alpha*A*x + y. * JX = KX IF (INCY.EQ.1) THEN DO 60 J = 1,N TEMP = ALPHA*X(JX) DO 50 I = 1,M Y(I) = Y(I) + TEMP*A(I,J) 50 CONTINUE JX = JX + INCX 60 CONTINUE ELSE DO 80 J = 1,N TEMP = ALPHA*X(JX) IY = KY DO 70 I = 1,M Y(IY) = Y(IY) + TEMP*A(I,J) IY = IY + INCY 70 CONTINUE JX = JX + INCX 80 CONTINUE END IF ELSE * * Form y := alpha*A**T*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 *> \brief \b DGEMM * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * * Definition: * =========== * * SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) * * .. Scalar Arguments .. * DOUBLE PRECISION ALPHA,BETA * INTEGER K,LDA,LDB,LDC,M,N * CHARACTER TRANSA,TRANSB * .. * .. Array Arguments .. * DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> 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**T, *> *> 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. *> \endverbatim * * Arguments: * ========== * *> \param[in] TRANSA *> \verbatim *> TRANSA is 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**T. *> *> TRANSA = 'C' or 'c', op( A ) = A**T. *> \endverbatim *> *> \param[in] TRANSB *> \verbatim *> TRANSB is 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**T. *> *> TRANSB = 'C' or 'c', op( B ) = B**T. *> \endverbatim *> *> \param[in] M *> \verbatim *> M is 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. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is 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. *> \endverbatim *> *> \param[in] K *> \verbatim *> K is 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. *> \endverbatim *> *> \param[in] ALPHA *> \verbatim *> ALPHA is DOUBLE PRECISION. *> On entry, ALPHA specifies the scalar alpha. *> \endverbatim *> *> \param[in] A *> \verbatim *> A is 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. *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is 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 ). *> \endverbatim *> *> \param[in] B *> \verbatim *> B is 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. *> \endverbatim *> *> \param[in] LDB *> \verbatim *> LDB is 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 ). *> \endverbatim *> *> \param[in] BETA *> \verbatim *> BETA is DOUBLE PRECISION. *> On entry, BETA specifies the scalar beta. When BETA is *> supplied as zero then C need not be set on input. *> \endverbatim *> *> \param[in,out] C *> \verbatim *> C is 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 ). *> \endverbatim *> *> \param[in] LDC *> \verbatim *> LDC is 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 ). *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2015 * *> \ingroup double_blas_level3 * *> \par Further Details: * ===================== *> *> \verbatim *> *> 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. *> \endverbatim *> * ===================================================================== SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) * * -- Reference BLAS level3 routine (version 3.6.0) -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2015 * * .. Scalar Arguments .. DOUBLE PRECISION ALPHA,BETA INTEGER K,LDA,LDB,LDC,M,N CHARACTER TRANSA,TRANSB * .. * .. Array Arguments .. DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) * .. * * ===================================================================== * * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB LOGICAL NOTA,NOTB * .. * .. Parameters .. DOUBLE PRECISION ONE,ZERO PARAMETER (ONE=1.0D+0,ZERO=0.0D+0) * .. * * 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 = LSAME(TRANSA,'N') NOTB = LSAME(TRANSB,'N') IF (NOTA) THEN NROWA = M NCOLA = K ELSE NROWA = K 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.LSAME(TRANSA,'C')) .AND. + (.NOT.LSAME(TRANSA,'T'))) THEN INFO = 1 ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND. + (.NOT.LSAME(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 TEMP = ALPHA*B(L,J) DO 70 I = 1,M C(I,J) = C(I,J) + TEMP*A(I,L) 70 CONTINUE 80 CONTINUE 90 CONTINUE ELSE * * Form C := alpha*A**T*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**T + 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 TEMP = ALPHA*B(J,L) DO 150 I = 1,M C(I,J) = C(I,J) + TEMP*A(I,L) 150 CONTINUE 160 CONTINUE 170 CONTINUE ELSE * * Form C := alpha*A**T*B**T + 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 *> \brief \b DISNAN tests input for NaN. * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download DISNAN + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * LOGICAL FUNCTION DISNAN( DIN ) * * .. Scalar Arguments .. * DOUBLE PRECISION DIN * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DISNAN returns .TRUE. if its argument is NaN, and .FALSE. *> otherwise. To be replaced by the Fortran 2003 intrinsic in the *> future. *> \endverbatim * * Arguments: * ========== * *> \param[in] DIN *> \verbatim *> DIN is DOUBLE PRECISION *> Input to test for NaN. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date September 2012 * *> \ingroup auxOTHERauxiliary * * ===================================================================== LOGICAL FUNCTION DISNAN( DIN ) * * -- LAPACK auxiliary routine (version 3.4.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * September 2012 * * .. Scalar Arguments .. DOUBLE PRECISION DIN * .. * * ===================================================================== * * .. External Functions .. LOGICAL DLAISNAN EXTERNAL DLAISNAN * .. * .. Executable Statements .. DISNAN = DLAISNAN(DIN,DIN) RETURN END logical function LSAME(A,B) character*1 A,B LSAME=A.eq.B return end *> \brief \b DSYR * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * * Definition: * =========== * * SUBROUTINE DSYR(UPLO,N,ALPHA,X,INCX,A,LDA) * * .. Scalar Arguments .. * DOUBLE PRECISION ALPHA * INTEGER INCX,LDA,N * CHARACTER UPLO * .. * .. Array Arguments .. * DOUBLE PRECISION A(LDA,*),X(*) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DSYR performs the symmetric rank 1 operation *> *> A := alpha*x*x**T + A, *> *> where alpha is a real scalar, x is an n element vector and A is an *> n by n symmetric matrix. *> \endverbatim * * Arguments: * ========== * *> \param[in] UPLO *> \verbatim *> UPLO is CHARACTER*1 *> On entry, UPLO specifies whether the upper or lower *> triangular part of the array A is to be referenced as *> follows: *> *> UPLO = 'U' or 'u' Only the upper triangular part of A *> is to be referenced. *> *> UPLO = 'L' or 'l' Only the lower triangular part of A *> is to be referenced. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> On entry, N specifies the order of the matrix A. *> N must be at least zero. *> \endverbatim *> *> \param[in] ALPHA *> \verbatim *> ALPHA is DOUBLE PRECISION. *> On entry, ALPHA specifies the scalar alpha. *> \endverbatim *> *> \param[in] X *> \verbatim *> X is 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. *> \endverbatim *> *> \param[in] INCX *> \verbatim *> INCX is INTEGER *> On entry, INCX specifies the increment for the elements of *> X. INCX must not be zero. *> \endverbatim *> *> \param[in,out] A *> \verbatim *> A is 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 part of the symmetric matrix and the strictly *> lower triangular part of A is not referenced. On exit, the *> upper triangular part of the array A is overwritten by the *> upper triangular part of the updated matrix. *> Before entry with UPLO = 'L' or 'l', the leading n by n *> lower triangular part of the array A must contain the lower *> triangular part of the symmetric matrix and the strictly *> upper triangular part of A is not referenced. On exit, the *> lower triangular part of the array A is overwritten by the *> lower triangular part of the updated matrix. *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is 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 ). *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2011 * *> \ingroup double_blas_level2 * *> \par Further Details: * ===================== *> *> \verbatim *> *> 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. *> \endverbatim *> * ===================================================================== SUBROUTINE DSYR(UPLO,N,ALPHA,X,INCX,A,LDA) * * -- Reference BLAS level2 routine (version 3.4.0) -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX,LDA,N CHARACTER UPLO * .. * .. Array Arguments .. DOUBLE PRECISION A(LDA,*),X(*) * .. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D+0) * .. * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I,INFO,IX,J,JX,KX * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * * Test the input parameters. * INFO = 0 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN INFO = 1 ELSE IF (N.LT.0) THEN INFO = 2 ELSE IF (INCX.EQ.0) THEN INFO = 5 ELSE IF (LDA.LT.MAX(1,N)) THEN INFO = 7 END IF IF (INFO.NE.0) THEN CALL XERBLA('DSYR ',INFO) RETURN END IF * * Quick return if possible. * IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN * * Set the start point in X if the increment is not unity. * 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 the triangular part * of A. * IF (LSAME(UPLO,'U')) THEN * * Form A when A is stored in upper triangle. * IF (INCX.EQ.1) THEN DO 20 J = 1,N IF (X(J).NE.ZERO) THEN TEMP = ALPHA*X(J) DO 10 I = 1,J A(I,J) = A(I,J) + X(I)*TEMP 10 CONTINUE END IF 20 CONTINUE ELSE JX = KX DO 40 J = 1,N IF (X(JX).NE.ZERO) THEN TEMP = ALPHA*X(JX) IX = KX DO 30 I = 1,J A(I,J) = A(I,J) + X(IX)*TEMP IX = IX + INCX 30 CONTINUE END IF JX = JX + INCX 40 CONTINUE END IF ELSE * * Form A when A is stored in lower triangle. * IF (INCX.EQ.1) THEN DO 60 J = 1,N IF (X(J).NE.ZERO) THEN TEMP = ALPHA*X(J) DO 50 I = J,N A(I,J) = A(I,J) + X(I)*TEMP 50 CONTINUE END IF 60 CONTINUE ELSE JX = KX DO 80 J = 1,N IF (X(JX).NE.ZERO) THEN TEMP = ALPHA*X(JX) IX = JX DO 70 I = J,N A(I,J) = A(I,J) + X(IX)*TEMP IX = IX + INCX 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF END IF * RETURN * * End of DSYR . * END *> \brief \b DLAISNAN tests input for NaN by comparing two arguments for inequality. * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * *> \htmlonly *> Download DLAISNAN + dependencies *> *> [TGZ] *> *> [ZIP] *> *> [TXT] *> \endhtmlonly * * Definition: * =========== * * LOGICAL FUNCTION DLAISNAN( DIN1, DIN2 ) * * .. Scalar Arguments .. * DOUBLE PRECISION DIN1, DIN2 * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> This routine is not for general use. It exists solely to avoid *> over-optimization in DISNAN. *> *> DLAISNAN checks for NaNs by comparing its two arguments for *> inequality. NaN is the only floating-point value where NaN != NaN *> returns .TRUE. To check for NaNs, pass the same variable as both *> arguments. *> *> A compiler must assume that the two arguments are *> not the same variable, and the test will not be optimized away. *> Interprocedural or whole-program optimization may delete this *> test. The ISNAN functions will be replaced by the correct *> Fortran 03 intrinsic once the intrinsic is widely available. *> \endverbatim * * Arguments: * ========== * *> \param[in] DIN1 *> \verbatim *> DIN1 is DOUBLE PRECISION *> \endverbatim *> *> \param[in] DIN2 *> \verbatim *> DIN2 is DOUBLE PRECISION *> Two numbers to compare for inequality. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date September 2012 * *> \ingroup auxOTHERauxiliary * * ===================================================================== LOGICAL FUNCTION DLAISNAN( DIN1, DIN2 ) * * -- LAPACK auxiliary routine (version 3.4.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * September 2012 * * .. Scalar Arguments .. DOUBLE PRECISION DIN1, DIN2 * .. * * ===================================================================== * * .. Executable Statements .. DLAISNAN = (DIN1.NE.DIN2) RETURN END *> \brief \b DGER * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * * Definition: * =========== * * SUBROUTINE DGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA) * * .. Scalar Arguments .. * DOUBLE PRECISION ALPHA * INTEGER INCX,INCY,LDA,M,N * .. * .. Array Arguments .. * DOUBLE PRECISION A(LDA,*),X(*),Y(*) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DGER performs the rank 1 operation *> *> A := alpha*x*y**T + 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. *> \endverbatim * * Arguments: * ========== * *> \param[in] M *> \verbatim *> M is INTEGER *> On entry, M specifies the number of rows of the matrix A. *> M must be at least zero. *> \endverbatim *> *> \param[in] N *> \verbatim *> N is INTEGER *> On entry, N specifies the number of columns of the matrix A. *> N must be at least zero. *> \endverbatim *> *> \param[in] ALPHA *> \verbatim *> ALPHA is DOUBLE PRECISION. *> On entry, ALPHA specifies the scalar alpha. *> \endverbatim *> *> \param[in] X *> \verbatim *> X is DOUBLE PRECISION array of dimension at least *> ( 1 + ( m - 1 )*abs( INCX ) ). *> Before entry, the incremented array X must contain the m *> element vector x. *> \endverbatim *> *> \param[in] INCX *> \verbatim *> INCX is INTEGER *> On entry, INCX specifies the increment for the elements of *> X. INCX must not be zero. *> \endverbatim *> *> \param[in] Y *> \verbatim *> Y is DOUBLE PRECISION array of dimension at least *> ( 1 + ( n - 1 )*abs( INCY ) ). *> Before entry, the incremented array Y must contain the n *> element vector y. *> \endverbatim *> *> \param[in] INCY *> \verbatim *> INCY is INTEGER *> On entry, INCY specifies the increment for the elements of *> Y. INCY must not be zero. *> \endverbatim *> *> \param[in,out] A *> \verbatim *> A is 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. On exit, A is *> overwritten by the updated matrix. *> \endverbatim *> *> \param[in] LDA *> \verbatim *> LDA is 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 ). *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2011 * *> \ingroup double_blas_level2 * *> \par Further Details: * ===================== *> *> \verbatim *> *> 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. *> \endverbatim *> * ===================================================================== SUBROUTINE DGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA) * * -- Reference BLAS level2 routine (version 3.4.0) -- * -- Reference BLAS is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX,INCY,LDA,M,N * .. * .. Array Arguments .. DOUBLE PRECISION A(LDA,*),X(*),Y(*) * .. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D+0) * .. * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I,INFO,IX,J,JY,KX * .. * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX * .. * * 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 readopt(tol,lred,rf) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) logical lex,lred character*80 s80 tol=1.0d-8 lred=.false. rf=1.0d-3 inquire(file='BUMATC6.opt',exist=lex) if(lex)then open(9,file='BUMATC6.opt') 1 read(9,80,end=99,err=99)s80 80 format(a80) c tol for B-matrix: if(s80(1:3).eq.'TOL')read(9,*)tol c tol for FF: if(s80(1:4).eq.'RFAC')read(9,*)rf c reduce Cartesian FF using rf: if(s80(1:4).eq.'REDU')read(9,*)lred goto 1 99 close(9) endif return end subroutine readff(MX,NAT3,FCAR) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) real*8 FCAR(MX,MX) WRITE(*,*)' FILE.FC ... INPUT FILE WITH' WRITE(*,*)' THE CARTESIAN FORCE FIELD' OPEN(20,FILE='FILE.FC',STATUS='OLD') N=NAT3 NTT= -4 NTE= 5 IF (NTE.GT.N) GO TO 120 110 NTT=NTT+5 DO 100 I=NTT,N NT=I IF (NT.GT.NTE) NT=NTE 100 READ(20,27) (FCAR(I,J),J=NTT,NT) 27 FORMAT(4X,5D14.6) NTE=NTE+5 IF (NTE.LE.N) GOTO 110 120 NTT=NTT+5 DO 130 I=NTT,N 130 READ(20,27) (FCAR(I,J),J=NTT,I) CLOSE(20) DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) DO 6 I=1,N DO 6 J=1,N C 6 FCAR(I,J)=FCAR(I,J)* 8.23D0 / 0.527D0 6 FCAR(I,J)=FCAR(I,J)*4.3597482D0/(0.529177249D0*0.529177249D0) c 1Hartree=4.3597482e-18 J; 1Bohr=0.529177249A write(6,*)'FILE.FC read' return end