* Auxiliary subroutines for MINL1, MININF, MINCF, MINCL1, MINCIN * Collected by H.B. Nielsen, IMM, DTU. 00.11.02 / 00.11.27 / 00.12.21 * * Jacobian checker: * SUBROUTINE CHECKD(FDF,N,M,X,STEPL,DIFF,INDX,A,A1,F,F1,F2,FAIL) ************************************************************************ * Check implementation of Jacobian of M functions of N variables * Hans Bruun Nielsen, IMM, DTU. 00.10.06 ************************************************************************ IMPLICIT NONE INTEGER N,M,INDX(6),FAIL, I,J DOUBLE PRECISION X(N),STEPL,DIFF(4),A(M,N),A1(M,N), & F(M),F1(M),F2(M), XJ,HB,HF,AF,AB,AE,ER EXTERNAL FDF INTRINSIC ABS, MAX C ... Initialize FAIL = 1 DO 10 I = 1, 4 10 DIFF(I) = 0D0 DO 20 I = 1, 3 20 INDX(I) = 0 CALL FDF(N,M,X,A,F) C ... Run through components of X DO 40 J = 1, N XJ = X(J) C ... Forward X(J) = XJ + STEPL HF = X(J) - XJ IF (HF .EQ. 0D0) RETURN CALL FDF(N,M,X,A1,F1) C ... Back X(J) = XJ - .5D0 * STEPL HB = X(J) - XJ IF (HB .EQ. 0D0) RETURN CALL FDF(N,M,X,A1,F2) C ... Run through components of f DO 30 I = 1, M DIFF(1) = MAX(DIFF(1), ABS(A(I,J))) AF = (F1(I) - F(I))/HF ER = AF - A(I,J) IF (ABS(ER) .GT. ABS(DIFF(2))) THEN DIFF(2) = ER INDX(1) = I INDX(2) = J ENDIF AB = (F2(I) - F(I))/HB ER = AB - A(I,J) IF (ABS(ER) .GT. ABS(DIFF(3))) THEN DIFF(3) = ER INDX(3) = I INDX(4) = J ENDIF C ... Extrapolated AE = (2D0*AB + AF)/3D0 ER = AE - A(I,J) IF (ABS(ER) .GT. ABS(DIFF(4))) THEN DIFF(4) = ER INDX(5) = I INDX(6) = J ENDIF 30 CONTINUE C ... Restore x(j) X(J) = XJ 40 CONTINUE FAIL = 0 RETURN END * * Simple output of a vector * SUBROUTINE PRVCTR(NAME,X,I1,I2,UNT) ************************************************************************ * Print on UNT elements I1 to I2 of X with name NAME * Hans Bruun Nielsen, Numerisk Institut, DTH. 89.09.28. ************************************************************************ CHARACTER*3 NAME INTEGER I1,I2,UNT, J,J1,J2 DOUBLE PRECISION X(*) C J2 = I1 - 1 10 IF (J2 .GE. I2) RETURN J1 = J2 + 1 J2 = J2 + 5 IF (J2 .GT. I2) J2 = I2 WRITE(UNT,'(1X,A,A,I3,A,I3,A,T18,1P1D10.3,1P4D12.3)') / NAME,'(',J1,'..',J2,') =',(X(J), J=J1,J2) GOTO 10 ************************** end of PRVCTR **************************** END * * LP subroutines etc. by J. Hald * SUBROUTINE L1LP(F,DF,C,DC,M,N,IC,LE,LI,XNMAX,XN,X,NACT, - KSET,ASET,U,R,DL,H,FUP,DLDF,CUP,DLDC,KSTATF,KSTATC,KSTATL,STEPF, - W,SEPS,ACCUM,FL1,ICONTR) C C THE SUBROUTINE SOLVES A LINEARLY CONSTRAINED C LINEAR L1 PROBLEM. C THE STARTING POINT MUST BE FEASIBLE. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION F(M),DF(M,N),C(IC),DC(IC,N),X(N),ASET(N),U(N,N), - R(N,N),DL(N),H(N),FUP(M),DLDF(M),CUP(IC),DLDC(IC),STEPF(M), - W(N) INTEGER KSET(N),KSTATF(M),KSTATC(IC),KSTATL(N) LOGICAL ACCUM,LINDEP C C INITIALIZE C LE1=LE+1 L=LE+LI LN=L+N LNN=LN+N XN=0D0 EPS=N*SEPS ACCUM=.FALSE. LINDEP=.FALSE. NCYC=0 ICONTR=0 FL1=1D73 DO 10 I=1,N KSTATL(I)=0 10 X(I)=0D0 DO 15 I=1,M KSTATF(I)=1 IF (F(I) .LT. 0D0) KSTATF(I)=-1 15 FUP(I)=F(I) C ACTIVATE INITIAL ACTIVE SET OF CONSTRAINTS NACTIN=NACT NACT=0 IF (L .EQ. 0) GOTO 100 DO 20 I=1,L KSTATC(I)=0 20 CUP(I)=C(I) IF (LE .EQ. 0) GOTO 50 DO 40 I=1,LE DO 30 J=1,N 30 R(J,I)=DC(I,J) KSET(I)=I 40 KSTATC(I)=1 CALL ADDCOL(U,R,N,NACT,LE,H,W,ACCUM,.FALSE.,EPS) IF (NACT .EQ. LE) GOTO 50 ICONTR=3 RETURN 50 IF (NACTIN .LT. LE1) GOTO 100 DO 70 K=LE1,NACTIN KK=KSET(K) IF (KK .LT. LE1 .OR. KK .GT. L) GOTO 70 NACT1=NACT+1 IF (NACT1 .GT. N) GOTO 100 DO 60 I=1,N 60 R(I,NACT1)=DC(KK,I) CALL UTTRNS(U,N,NACT,ACCUM,R(1,NACT1),W) CALL ADDCOL(U,R,N,NACT,1,H,W,ACCUM,.FALSE.,EPS) IF (NACT .LT. NACT1) GOTO 70 EPS=EPS+2D0*SEPS KSET(NACT1)=KK KSTATC(KK)=1 70 CONTINUE C CALCULATE THE RESTRICTED GRADIENT. 100 FL1=0D0 DO 105 I=1,N 105 H(I)=0D0 DO 140 I=1,M FL1=FL1+DABS(FUP(I)) K=KSTATF(I) IF (K .GT. -1) GOTO 120 DO 110 J=1,N 110 H(J)=H(J)-DF(I,J) 120 IF (K .LT. 1) GOTO 140 DO 130 J=1,N 130 H(J)=H(J)+DF(I,J) 140 CONTINUE C TRANSFORM THE RESTRICTED GRADIENT. DO 150 I=1,N 150 DL(I)=-H(I) CALL UTTRNS(U,N,NACT,ACCUM,DL,W) C CHECK IF THE CURRENT POINT IS A DEAD POINT. 200 T=0D0 NACT1=NACT+1 DLN2=0D0 IF (NACT .EQ. N) GOTO 220 DO 210 I=NACT1,N 210 DLN2=DLN2+DL(I)**2 220 IF (NACT1 .EQ. 1) GOTO 240 DO 230 I=1,NACT 230 T=T+DL(I)**2 240 T=T+DLN2 IF (T .EQ. 0D0) GOTO 900 IF (DLN2 .GT. EPS*EPS*T) GOTO 395 C DEAD POINT - CALCULATE MULTIPLIERS. IF (NACT .EQ. LE) GOTO 900 CALL RSOLV(R,N,NACT,DL,ASET) AMAX=-1D73 DO 250 K=LE1,NACT A=ASET(K) IF (KSET(K) .GT. LNN) A=DABS(A)-1D0 IF (A .LE. AMAX) GOTO 250 KMAX=K AMAX=A 250 CONTINUE IF (AMAX .LE. 0D0) GOTO 900 C DELETE A FUNCTION FROM THE ACTIVE SET. KK=KSET(KMAX) IF (KK .LE. LNN) GOTO 310 KK=KK-LNN A=ASET(KMAX) IF (A .LT. 0D0) GOTO 280 KSTATF(KK)=1 DO 260 I=1,N 260 H(I)=H(I)+DF(KK,I) DO 270 I=1,KMAX 270 DL(I)=DL(I)-R(I,KMAX) GOTO 330 280 KSTATF(KK)=-1 DO 290 I=1,N 290 H(I)=H(I)-DF(KK,I) DO 300 I=1,KMAX 300 DL(I)=DL(I)+R(I,KMAX) GOTO 330 C DEACTIVATE A CONSTRAINT. 310 IF (KK .GT. L) GOTO 320 KSTATC(KK)=0 GOTO 330 320 IF (KK .GT. LN) GOTO 325 KSTATL(KK-L)=0 GOTO 330 325 KSTATL(KK-LN)=0 C REMOVE THE CORRESPONDING COLUMN IN R. 330 IF (ACCUM) GOTO 340 ACCUM=.TRUE. CALL HACCUM(U,N,NACT,W) 340 CALL DELCOL(KMAX,U,R,N,NACT,DL,.TRUE.) EPS=EPS+2D0*SEPS IF (KMAX .GT. NACT) GOTO 360 DO 350 I=KMAX,NACT 350 KSET(I)=KSET(I+1) C REMOVE LINEAR DEPENDENCE LABELS. 360 IF (.NOT. LINDEP) GOTO 200 LINDEP=.FALSE. DO 370 I=1,M IF (KSTATF(I) .EQ. 2) KSTATF(I)=1 370 IF (KSTATF(I) .EQ. -2) KSTATF(I)=-1 DO 380 I=1,N 380 IF (KSTATL(I) .EQ. 2) KSTATL(I)=0 IF (LI .EQ. 0) GOTO 200 DO 390 I=LE1,L 390 IF (KSTATC(I) .EQ. 2) KSTATC(I)=0 GOTO 200 C CALCULATE THE DIRECTION OF SEARCH. 395 IF (NACT .EQ. 0) GOTO 405 DO 400 I=1,NACT 400 DL(I)=0D0 405 CALL UTRNS(U,N,NACT,ACCUM,DL,W) DLH=0D0 DO 410 I=1,N 410 DLH=DLH+DL(I)*H(I) C PROJECT GRADIENTS ON THE DIRECTION OF SEARCH. CALL MATVEC(DF,DL,M,M,N,DLDF) CALL MATVEC(DC,DL,IC,L,N,DLDC) C CALCULATE THE STEP LENGTHS. DO 420 I=1,M S=DLDF(I) K=KSTATF(I) STEPF(I)=1D73 IF (IABS(K) .NE. 1 .OR. K*S .GE. -1D-50) GOTO 420 STEPF(I)=-FUP(I)/S 420 CONTINUE 425 SLP=1D73 SLN=1D73 DO 440 I=1,N IF (KSTATL(I) .NE. 0) GOTO 440 S=DL(I) IF (S*S .LT. EPS*EPS*DLN2) GOTO 440 T=(XNMAX-X(I))/S IF (S .LT. 0D0 .OR. T .GE. SLP) GOTO 430 SLP=T KLP=I 430 T=-(XNMAX+X(I))/S IF (S .GT. 0D0 .OR. T .GE. SLN) GOTO 440 SLN=T KLN=I 440 CONTINUE SC=1D73 IF (LI .EQ. 0) GOTO 460 DO 450 I=LE1,L S=DLDC(I) IF (KSTATC(I) .NE. 0 .OR. S .GE. -1D-50) GOTO 450 T=-CUP(I)/S IF (T .GE. SC) GOTO 450 KC=I SC=T 450 CONTINUE 460 STEPLC=DMIN1(SLP,SLN) STEPLC=DMIN1(STEPLC,SC) NACT1=NACT+1 C DO LINE SEARCH. 500 SF=1D73 DO 510 I=1,M S=STEPF(I) IF (S .GE. SF) GOTO 510 KF=I SF=S 510 CONTINUE STEP=DMIN1(STEPLC,SF) IF (SF .GT. STEP) GOTO 630 T=DLH-2*KSTATF(KF)*DLDF(KF) IF (T .GE. 0D0) GOTO 600 DLH=T KSTATF(KF)=-KSTATF(KF) STEPF(KF)=1D73 GOTO 500 C ACTIVATE A FUNCTION. 600 DO 610 I=1,N 610 R(I,NACT1)=DF(KF,I) CALL UTTRNS(U,N,NACT,ACCUM,R(1,NACT1),W) CALL ADDCOL(U,R,N,NACT,1,DL,W,ACCUM,.FALSE.,EPS) IF (NACT .EQ. NACT1) GOTO 620 KSTATF(KF)=2*KSTATF(KF) STEPF(KF)=1D73 LINDEP=.TRUE. GOTO 500 620 KSET(NACT)=KF+LNN KSTATF(KF)=0 GOTO 700 C - OR A CONSTRAINT. 630 IF (SC .GT. STEP) GOTO 660 DO 640 I=1,N 640 R(I,NACT1)=DC(KC,I) CALL UTTRNS(U,N,NACT,ACCUM,R(1,NACT1),W) CALL ADDCOL(U,R,N,NACT,1,DL,W,ACCUM,.FALSE.,EPS) IF (NACT .EQ. NACT1) GOTO 650 KSTATC(KC)=2 LINDEP=.TRUE. GOTO 425 650 KSTATC(KC)=1 KSET(NACT)=KC GOTO 700 C - OR A BOUND. 660 DO 670 I=1,N 670 R(I,NACT1)=0D0 IF (SLN .EQ. SLP) SLN=2D73 IF (STEP .EQ. SLN) R(KLN,NACT1)=1D0 IF (STEP .EQ. SLP) R(KLP,NACT1)=-1D0 CALL UTTRNS(U,N,NACT,ACCUM,R(1,NACT1),W) CALL ADDCOL(U,R,N,NACT,1,DL,W,ACCUM,.FALSE.,EPS) IF (NACT .EQ. NACT1) GOTO 680 LINDEP=.TRUE. IF (STEP .EQ. SLN) KSTATL(KLN)=2 IF (STEP .EQ. SLP) KSTATL(KLP)=2 GOTO 425 680 IF (STEP .EQ. SLN) GOTO 690 KSET(NACT)=KLP+L KSTATL(KLP)=1 GOTO 700 690 KSET(NACT)=KLN+LN KSTATL(KLN)=-1 700 EPS=EPS+2D0*SEPS C INTRODUCE THE NEW POINT. XN=0D0 XXN=0D0 NCYC=NCYC+1 IF (NCYC .GT. 2*N) GOTO 900 DO 800 I=1,N S=STEP*DL(I) XXN=DMAX1(XXN,DABS(S)) X(I)=X(I)+S 800 XN=DMAX1(XN,DABS(X(I))) IF (XXN .GT. EPS*XN) NCYC=0 DO 810 I=1,M 810 FUP(I)=FUP(I)+STEP*DLDF(I) IF (L .EQ. 0) GOTO 100 DO 820 I=1,L 820 CUP(I)=CUP(I)+STEP*DLDC(I) GOTO 100 C RETURN 900 IF (NACT .LT. N) ICONTR=2 DO 910 I=1,NACT 910 IF (KSET(I) .GT. L .AND. KSET(I) .LE. LNN) ICONTR=1 IF (ICONTR .EQ. 1) XN=XNMAX IF (NCYC .GT. 2*N) ICONTR=4 RETURN END *** SUBROUTINE FEASI(C,DC,IC,LE,LI,N,X,NACT,KSET,ASET,U, - R,DL,RIGHT,CUP,DLDC,W,KSTAT,ICONTR,ACCUM,SEPS) C C THE SUBROUTINE FINDS A FEASIBLE POINT FOR A SET OF LINEAR C EQUALITY AND INEQUALITY CONSTRAINTS. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION C(IC),DC(IC,N),X(N),ASET(N),U(N,N),R(N,N),DL(N), - RIGHT(N),CUP(IC),DLDC(IC),W(N) INTEGER KSET(N),KSTAT(IC) LOGICAL ACCUM,OBJECT C INITIALIZE EPS=(N+10)*SEPS ACCUM=.FALSE. NACTIN=NACT NACT=0 LE1=LE+1 LELI=LE+LI DO 10 I=1,N 10 X(I)=0D0 ICONTR=0 IF (LELI .EQ. 0) RETURN DO 15 I=1,LELI 15 KSTAT(I)=0 C MAKE ACTIVE THE EQUALITY CONSTRAINTS PLUS OTHER CONSTRAINTS C AS DEFINED IN KSET. IF (LE .EQ. 0) GOTO 35 IF (LE .GT. N) GOTO 600 DO 30 I=1,LE RIGHT(I)=-C(I) DO 20 J=1,N 20 R(J,I)=DC(I,J) KSET(I)=I 30 KSTAT(I)=1 CALL ADDCOL(U,R,N,NACT,LE,RIGHT,W,ACCUM,.FALSE.,EPS) IF (NACT .LT. LE) GOTO 600 35 IF (NACTIN .LT. 1) GOTO 60 DO 50 K=1,NACTIN KK=KSET(K) IF(KK .LT. LE1 .OR. KK .GT. LELI) GOTO 50 NACT1=NACT+1 IF (NACT1 .GT. N) GOTO 60 DO 40 I=1,N 40 R(I,NACT1)=DC(KK,I) CALL UTTRNS(U,N,NACT,ACCUM,R(1,NACT1),W) CALL ADDCOL(U,R,N,NACT,1,RIGHT,W,ACCUM,.FALSE.,EPS) IF (NACT .LT. NACT1) GOTO 50 RIGHT(NACT1)=-C(KK) KSET(NACT1)=KK KSTAT(KK)=1 50 CONTINUE 60 CALL RTSOLV(R,N,NACT,RIGHT,X) IF (NACT .EQ. N) GOTO 80 NACT1=NACT+1 DO 70 I=NACT1,N 70 X(I)=0D0 80 CALL UTRNS(U,N,NACT,ACCUM,X,W) C UPDATE THE CONSTRAINTS. CALL MATVEC(DC,X,IC,LELI,N,CUP) DO 100 I=1,LELI 100 CUP(I)=CUP(I)+C(I) C INITIALIZE INEQUALITY CONSTRAINT LOOP. DO 120 I=LE1,LELI 120 IF( CUP(I) .LT. 0D0 .AND. KSTAT(I) .EQ. 0) KSTAT(I)=-1 C ACTIVATE VIOLATED INEQUALITY CONSTRAINTS ONE BY ONE. C USE THE STRONGEST VIOLATED AS OBJECTIVE CONSTRAINT. 150 FMIN=1D73 DO 160 I=LE1,LELI IF (KSTAT(I) .NE. -1) GOTO 160 IF (CUP(I) .GE. FMIN) GOTO 160 FMIN=CUP(I) NEW=I 160 CONTINUE IF (FMIN .GE. 0D0) RETURN DO 170 I=1,N 170 RIGHT(I)=DC(NEW,I) CALL UTTRNS(U,N,NACT,ACCUM,RIGHT,W) KSTAT(NEW)=1 C DEAD POINT? 200 IF (NACT .EQ. N) GOTO 240 T=0D0 DLN2=0D0 IF (NACT .EQ. 0) GOTO 220 DO 210 I=1,NACT T=T+RIGHT(I)**2 210 DL(I)=0D0 220 NACT1=NACT+1 DO 230 I=NACT1,N DLN2=DLN2+RIGHT(I)**2 230 DL(I)=RIGHT(I) T=T+DLN2 IF (T.GT. 1D-40 .AND. DLN2 .GT. EPS*EPS*T) GOTO 300 C - YES - CALCULATE MULTIPLIERS AND SKIP AN ACTIVE C CONSTRAINT IF POSSIBLE. 240 IF (NACT .EQ. LE) GOTO 600 CALL RSOLV(R,N,NACT,RIGHT,ASET) AMAX=-1D73 DO 250 I=LE1,NACT IF (ASET(I) .LT. AMAX) GOTO 250 K=I AMAX=ASET(I) 250 CONTINUE IF (AMAX .LE. 0D0) GOTO 600 KSTAT(KSET(K))=0 DO 260 I=LE1,LELI 260 IF (KSTAT(I) .EQ. -2) KSTAT(I)=0 IF (ACCUM) GOTO 270 ACCUM=.TRUE. CALL HACCUM(U,N,NACT,W) 270 CALL DELCOL(K,U,R,N,NACT,RIGHT,.TRUE.) IF (K .GT. NACT) GOTO 200 DO 280 I=K,NACT 280 KSET(I)=KSET(I+1) GOTO 200 C - NO - CALCULATE STEP DIRECTION. 300 CALL UTRNS(U,N,NACT,ACCUM,DL,W) C PROJECT GRADIENTS ON STEP DIRECTION. CALL MATVEC(DC,DL,IC,LELI,N,DLDC) C CALCULATE STEP LENGTH ANES TO MAKE THE OBJECTIVE CONSTRAINT C EQUAL ZERO, AND CALCULATE THE STEP LENGTH AMIN TO THE C NEAREST INACTIVE CONSTRAINT UNDER CONSIDERATION. ANES=-CUP(NEW)/DLN2 310 AMIN=1D73 DO 320 I=LE1,LELI IF (KSTAT(I) .NE. 0) GOTO 320 T=DLDC(I) IF (T .GE. -1D-50) GOTO 320 T=-CUP(I)/T IF (T .GT. AMIN) GOTO 320 AMIN=T K=I 320 CONTINUE C WILL THE OBJECTIVE CONSTRAINT GET ACTIVE C IF NOT, MAKE ACTIVE THE CLOSEST. OBJECT= ANES .LE. AMIN ALFA=DMIN1(AMIN,ANES) NACT1=NACT+1 IF (OBJECT) GOTO 370 DO 350 I=1,N 350 R(I,NACT1)=DC(K,I) CALL UTTRNS(U,N,NACT,ACCUM,R(1,NACT1),W) CALL ADDCOL(U,R,N,NACT,1,RIGHT,W,ACCUM,.TRUE.,EPS) IF (NACT1 .EQ. NACT) GOTO 360 KSTAT(K)=-2 GOTO 310 360 KSTAT(K)=1 KSET(NACT)=K C TAKE THE STEP. 370 IF (ALFA .EQ. 0D0) GOTO 410 DO 380 I=1,N 380 X(I)=X(I)+ALFA*DL(I) DO 400 I=1,LELI T=CUP(I)+ALFA*DLDC(I) IF (KSTAT(I) .EQ. -1 .AND. T .GE. 0D0) KSTAT(I)=0 400 CUP(I)=T 410 IF (.NOT. OBJECT) GOTO 200 C ACTIVATE THE OBJECTIVE CONSTRAINT. DO 500 I=1,N 500 R(I,NACT1)=RIGHT(I) CALL ADDCOL(U,R,N,NACT,1,RIGHT,W,ACCUM,.FALSE.,EPS) IF (NACT .EQ. NACT1) GOTO 510 KSTAT(NEW)=0 GOTO 150 510 KSET(NACT)=NEW GOTO 150 C NO FEASIBLE POINTS. 600 ICONTR=3 RETURN END *** SUBROUTINE LINSYS(A,B,IDIM,N,NR,SEPS) C THE SUBROUTINE SOLVES A SYSTEM OF LINEAR EQUATIONS C USING GAUSSIAN ELIMINATION IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION A(IDIM,IDIM),B(N) NR=0 C A IS CONSIDERED TO BE OF RANK K-1 IF THE ABSOLUTE VALUE C OF THE K'TH PIVOT IS LESS THAN K*SEPS. IF (N-1) 200,1,10 1 IF (DABS(A(1,1)) .LT. 1D-50) RETURN NR=1 B(1)=B(1)/A(1,1) RETURN C EQUILIBRATION IN THE INFINITY NORM. 10 DO 30 I=1,N AM=DABS(A(I,1)) DO 20 J=2,N S=DABS(A(I,J)) 20 IF (AM .LT. S) AM=S IF (AM .LT. 1D-50) AM=1D0 B(I)=B(I)/AM DO 30 J=1,N 30 A(I,J)=A(I,J)/AM C ELIMINATION N1=N-1 DO 80 K=1,N1 NR=K-1 C FIND PIVOTAL ROW. AM=DABS(A(K,K)) I0=K K1=K+1 DO 40 I=K1,N S=DABS(A(I,K)) IF (S .LE. AM) GOTO 40 AM=S I0=I 40 CONTINUE IF (AM .LT. 2*K*SEPS) RETURN IF (I0 .EQ. K) GOTO 60 C INTERCHANGE EQUATIONS K AND I0 DO 50 J=K,N S=A(K,J) A(K,J)=A(I0,J) 50 A(I0,J)=S S=B(K) B(K)=B(I0) B(I0)=S C STORE PIVOT IN AM AND ELIMINATE IN ROWS K+1 TO N. 60 AM=A(K,K) DO 80 I=K1,N S=A(I,K)/AM DO 70 J=K1,N 70 A(I,J)=A(I,J)-S*A(K,J) 80 B(I)=B(I)-S*B(K) NR=N1 IF (DABS(A(N,N)) .LT. 2*N*SEPS) RETURN C A HAS FULL RANK NR=N C BACK SUBSTITUTION. B(N)=B(N)/A(N,N) K=N DO 100 I=2,N K1=K K=K-1 S=B(K) DO 90 J=K1,N 90 S=S-A(K,J)*B(J) 100 B(K)=S/A(K,K) 200 RETURN END *** SUBROUTINE BFGS(B,N,Y,XX,W,SEPS) C UPDATES A HESSIAN APPROXIMATION USING BFGS-FORMULA. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION B(N,N),Y(N),XX(N),W(N) EPS=(N+10)*SEPS DO 20 I=1,N T=0D0 DO 10 J=1,N 10 T=T+B(I,J)*XX(J) 20 W(I)=T YXX=0D0 WXX=0D0 DO 30 I=1,N YXX=YXX+Y(I)*XX(I) 30 WXX=WXX+W(I)*XX(I) IF (WXX .LT. 1D-40 .OR. YXX .LT. EPS*WXX) RETURN DO 40 I=1,N 40 B(I,I)=B(I,I)+Y(I)**2/YXX-W(I)**2/WXX IF(N.EQ.1) RETURN DO 50 I=2,N I1=I-1 DO 50 J=1,I1 B(I,J)=B(I,J)+Y(I)*Y(J)/YXX-W(I)*W(J)/WXX 50 B(J,I)=B(I,J) RETURN END *** SUBROUTINE ADDCOL(U,R,N,KCOL,KNEW,RIGHT,W,ACCUM,LRIGHT,EPS) C UPDATES HOUSEHOLDER FACTORIZATION. C THE NEW COLUMNS MUST HAVE BEEN TRANSFORMED C AS RIGHTHAND SIDES. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(N,N),R(N,N),RIGHT(N),W(N) LOGICAL ACCUM,LRIGHT K1=KCOL+1 K2=KCOL+KNEW C COLUMN LOOP STARTS HERE DO 200 K=K1,K2 S=0D0 T=0D0 IF (K .EQ. 1) GOTO 20 KK=K-1 DO 10 I=1,KK 10 T=T+R(I,K)**2 20 DO 30 I=K,N 30 S=S+R(I,K)**2 T=T+S T=DSQRT(T) S=DSQRT(S) C RETURN IF THE NEW COLUMN DEPENDS LINEARLY ON THE C PRECEEDING COLUMNS IF (T .EQ. 0D0) RETURN IF (S .LT. T*EPS) RETURN C PERFORM HOUSEHOLDER TRANSFORMATION TT=R(K,K) T=DABS(TT) ALFA=DSQRT(S*(S+T)) BETA=-DSIGN(S,TT) R(K,K)=BETA W(K)=(TT-BETA)/ALFA IF (K .EQ. N) GOTO 80 KK=K+1 DO 40 I=KK,N 40 W(I)=R(I,K)/ALFA C TRANSFORM THE REMAINING COLUMNS IF (K .EQ. K2) GOTO 80 DO 70 J=KK,K2 T=0D0 DO 50 I=K,N 50 T=T+W(I)*R(I,J) DO 60 I=K,N 60 R(I,J)=R(I,J)-T*W(I) 70 CONTINUE C TRANSFORM THE RIGHTHAND SIDE 80 IF (.NOT. LRIGHT) GOTO 110 T=0D0 DO 90 I=K,N 90 T=T+W(I)*RIGHT(I) DO 100 I=K,N 100 RIGHT(I)=RIGHT(I)-T*W(I) C ACCUMULATE THE TRANSFORMATIONS IN U C U MUST HAVE BEEN INITIALIZED 110 IF (ACCUM) GOTO 130 DO 120 I=K,N 120 U(I,K)=W(I) GOTO 200 130 DO 160 I=1,N T=0D0 DO 140 J=K,N 140 T=T+U(I,J)*W(J) DO 150 J=K,N 150 U(I,J)=U(I,J)-T*W(J) 160 CONTINUE 200 KCOL=KCOL+1 RETURN END *** SUBROUTINE DELCOL(K,U,R,N,KCOL,RIGHT,LRIGHT) C DELETES COLUMN NO. K IN THE FACTORIZED MATRIX. C K MUST SATISFY 1.LE.K.LE.KCOL C U MUST HAVE BEEN ACCUMULATED. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(N,N),R(N,N),RIGHT(N) LOGICAL LRIGHT C DELETE COLUMN NUMBER K KCOL=KCOL-1 IF (K .GT. KCOL) RETURN DO 10 J=K,KCOL J1=J+1 DO 10 I=1,J1 10 R(I,J)=R(I,J1) C TRANSFORM TO UPPER TRIANGULAR FORM C USING STANDARD GIVENS TRANSFORMATIONS DO 100 KK=K,KCOL K1=KK+1 X=R(KK,KK) Y=R(K1,KK) A=DSQRT(X*X+Y*Y) C=X/A S=Y/A R(KK,KK)=C*X+S*Y IF (KK .EQ. KCOL) GOTO 60 DO 50 J=K1,KCOL X=R(KK,J) Y=R(K1,J) R(KK,J)=C*X+S*Y 50 R(K1,J)=C*Y-S*X 60 IF (.NOT. LRIGHT) GOTO 70 X=RIGHT(KK) Y=RIGHT(K1) RIGHT(KK)=C*X+S*Y RIGHT(K1)=C*Y-S*X C ACCUMULATE THE TRANSFORMATIONS 70 DO 80 I=1,N X=U(I,KK) Y=U(I,K1) U(I,KK)=C*X+S*Y 80 U(I,K1)=C*Y-S*X 100 CONTINUE RETURN END *** SUBROUTINE UTTRNS(U,N,KCOL,ACCUM,R,W) C TRANSFORM THE VECTOR R AS A RIGHTHAND SIDE. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(N,N),R(N),W(N) LOGICAL ACCUM C IF THE TRANSFORMATIONS HAVE BEEN ACCUMULATED C DO SIMPLE MATRIX-MULTIPLICATION. C ELSE TRANSFORM LIKE RIGHTHAND SIDES. IF (ACCUM) GOTO 100 IF (KCOL .EQ. 0) RETURN DO 30 K=1,KCOL T=0D0 DO 10 I=K,N 10 T=T+R(I)*U(I,K) DO 20 I=K,N 20 R(I)=R(I)-T*U(I,K) 30 CONTINUE RETURN 100 DO 110 I=1,N 110 W(I)=R(I) DO 130 K=1,N T=0D0 DO 120 I=1,N 120 T=T+U(I,K)*W(I) 130 R(K)=T RETURN END *** SUBROUTINE UTRNS(U,N,KCOL,ACCUM,R,W) C TRANSFORM THE VECTOR R OPPOSITE A RIGHTHAND SIDE. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(N,N),R(N),W(N) LOGICAL ACCUM K1=KCOL+1 IF (ACCUM) GOTO 100 IF (KCOL .EQ. 0) RETURN DO 30 KK=1,KCOL K=K1-KK T=0D0 DO 10 J=K,N 10 T=T+U(J,K)*R(J) DO 20 J=K,N 20 R(J)=R(J)-T*U(J,K) 30 CONTINUE RETURN 100 DO 110 I=1,N 110 W(I)=R(I) DO 130 I=1,N T=0D0 DO 120 J=1,N 120 T=T+U(I,J)*W(J) 130 R(I)=T RETURN END *** SUBROUTINE RSOLV(R,N,KCOL,RIGHT,X) C PERFORM BACK SUBSTITUTION ON RIGHT. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION R(N,N),RIGHT(N),X(N) C CALCULATE ALFA USING BACK SUBSTITUTION ON R K=KCOL K1=K+1 10 IF (K .EQ. 0) RETURN T=RIGHT(K) IF (K1 .GT. KCOL) GOTO 30 DO 20 J=K1,KCOL 20 T=T-X(J)*R(K,J) 30 X(K)=T/R(K,K) K1=K K=K-1 GOTO 10 END *** SUBROUTINE RTSOLV(R,N,KCOL,RIGHT,X) C PERFORM BACK SUBSTITUTION ON RIGHT USING THE C TRANSPOSED TRIANGULAR MATRIX. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION R(N,N),RIGHT(N),X(N) IF (KCOL .EQ. 0) RETURN X(1)=RIGHT(1)/R(1,1) IF (KCOL .EQ. 1) RETURN DO 20 I=2,KCOL I1=I-1 T=RIGHT(I) DO 10 J=1,I1 10 T=T-X(J)*R(J,I) 20 X(I)=T/R(I,I) RETURN END *** SUBROUTINE HACCUM(U,N,KCOL,W) C ACCUMULATES HOUSEHOLDER VECORS STORED IN LOWER TRIANGLE C OF THE FIRST KCOL COLUMNS OF U IN AN ORTHONORMAL MATRIX U. C THE HOUSEHOLDER VECTORS MUST HAVE TWO NORM EQUAL TO TWO. C KCOL.GE.1 . IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION U(N,N),W(N) C INITIALIZE USING LAST TRANSFORMATION K1=KCOL+1 DO 10 I=KCOL,N 10 W(I)=U(I,KCOL) DO 20 I=KCOL,N 20 U(I,I)=1D0-W(I)**2 IF (KCOL .EQ. N) GOTO 40 DO 30 I=K1,N I1=I-1 T=W(I) DO 30 J=KCOL,I1 S=-T*W(J) U(I,J)=S 30 U(J,I)=S 40 IF (KCOL .EQ. 1) RETURN C ACCUMULATE REMAINING TRANSFORMATIONS DO 100 KK=2,KCOL K=K1-KK DO 50 I=K,N 50 W(I)=U(I,K) T=W(K) KP1=K+1 U(K,K)=1D0-T**2 DO 60 I=KP1,N 60 U(I,K)=-T*W(I) DO 90 L=KP1,N S=0D0 DO 70 I=KP1,N 70 S=S+W(I)*U(I,L) U(K,L)=-T*S DO 80 I=KP1,N 80 U(I,L)=U(I,L)-S*W(I) 90 CONTINUE 100 CONTINUE RETURN END *** SUBROUTINE LIMIT(XNMAX2,X,XN2,P,PN2,ALFA,N) C LIMIT THE STEP LENGTH ALFA. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N),P(N) XTP=0D0 DO 10 I=1,N 10 XTP=XTP+X(I)*P(I) B=XTP/PN2 T=DSQRT(B*B+(XNMAX2-XN2)/PN2) AP=T-B AM=-T-B IF (ALFA .GT. AP) ALFA=AP IF (ALFA .LT. AM) ALFA=AM RETURN END *** SUBROUTINE MATVEC(A,X,IA,M,N,Y) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(IA,N),X(N),Y(IA) IF (M .EQ. 0) RETURN DO 20 I=1,M T=0D0 DO 10 J=1,N 10 T=T+A(I,J)*X(J) 20 Y(I)=T RETURN END C ************************************************************************ * Selection BLAS subroutines and functions, * Level 0 : LSAME, XERBLA * Level 1 : DAXPY, DCOPY, DDOT, DNRM2, DSCAL, IDAMAX * Level 2 : DGEMV, DTRSV * Level 3 : DGEMM, DSYR * Below we give the version obtained from * http://www.netlib.org/blas/blas.tgz ************************************************************************ LOGICAL FUNCTION LSAME( CA, CB ) * * -- LAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * January 31, 1994 * * .. Scalar Arguments .. CHARACTER CA, CB * .. * * Purpose * ======= * * LSAME returns .TRUE. if CA is the same letter as CB regardless of * case. * * Arguments * ========= * * CA (input) CHARACTER*1 * CB (input) CHARACTER*1 * CA and CB specify the single characters to be compared. * * ===================================================================== * * .. Intrinsic Functions .. INTRINSIC ICHAR * .. * .. Local Scalars .. INTEGER INTA, INTB, ZCODE * .. * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA.EQ.CB IF( LSAME ) $ RETURN * * Now test for equivalence if both characters are alphabetic. * ZCODE = ICHAR( 'Z' ) * * Use 'Z' rather than 'A' so that ASCII can be detected on Prime * machines, on which ICHAR returns a value with bit 8 set. * ICHAR('A') on Prime machines returns 193 which is the same as * ICHAR('A') on an EBCDIC machine. * INTA = ICHAR( CA ) INTB = ICHAR( CB ) * IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN * * ASCII is assumed - ZCODE is the ASCII code of either lower or * upper case 'Z'. * IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 * ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN * * EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or * upper case 'Z'. * IF( INTA.GE.129 .AND. INTA.LE.137 .OR. $ INTA.GE.145 .AND. INTA.LE.153 .OR. $ INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR. $ INTB.GE.145 .AND. INTB.LE.153 .OR. $ INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 * ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN * * ASCII is assumed, on Prime machines - ZCODE is the ASCII code * plus 128 of either lower or upper case 'Z'. * IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LSAME = INTA.EQ.INTB * * RETURN * * End of LSAME * END ************************************************************************ SUBROUTINE XERBLA( SRNAME, INFO ) * * -- LAPACK auxiliary routine (preliminary version) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * February 29, 1992 * * .. Scalar Arguments .. CHARACTER*6 SRNAME INTEGER INFO * .. * * Purpose * ======= * * XERBLA is an error handler for the LAPACK routines. * It is called by an LAPACK routine if an input parameter has an * invalid value. A message is printed and execution stops. * * Installers may consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Arguments * ========= * * SRNAME (input) CHARACTER*6 * The name of the routine which called XERBLA. * * INFO (input) INTEGER * The position of the invalid parameter in the parameter list * of the calling routine. * * WRITE( *, FMT = 9999 )SRNAME, INFO * STOP * 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ', $ 'an illegal value' ) * * End of XERBLA * END ************************************************************************ * Level 1 ************************************************************************ subroutine daxpy(n,da,dx,incx,dy,incy) c c constant times a vector plus a vector. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),da integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dy(i) + da*dx(i) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 dy(i) = dy(i) + da*dx(i) dy(i + 1) = dy(i + 1) + da*dx(i + 1) dy(i + 2) = dy(i + 2) + da*dx(i + 2) dy(i + 3) = dy(i + 3) + da*dx(i + 3) 50 continue return end ************************************************************************ subroutine dcopy(n,dx,incx,dy,incy) c c copies a vector, x, to a vector, y. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*) integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dy(iy) = dx(ix) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,7) if( m .eq. 0 ) go to 40 do 30 i = 1,m dy(i) = dx(i) 30 continue if( n .lt. 7 ) return 40 mp1 = m + 1 do 50 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) 50 continue return end ************************************************************************ double precision function ddot(n,dx,incx,dy,incy) c c forms the dot product of two vectors. c uses unrolled loops for increments equal to one. c jack dongarra, linpack, 3/11/78. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dy(*),dtemp integer i,incx,incy,ix,iy,m,mp1,n c ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + * dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end ************************************************************************ DOUBLE PRECISION FUNCTION DNRM2 ( N, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, N * .. Array Arguments .. DOUBLE PRECISION X( * ) * .. * * DNRM2 returns the euclidean norm of a vector via the function * name, so that * * DNRM2 := sqrt( x'*x ) * * * * -- This version written on 25-October-1982. * Modified on 14-October-1993 to inline the call to DLASSQ. * Sven Hammarling, Nag Ltd. * * * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. INTEGER IX DOUBLE PRECISION ABSXI, NORM, SCALE, SSQ * .. Intrinsic Functions .. INTRINSIC ABS, SQRT * .. * .. Executable Statements .. IF( N.LT.1 .OR. INCX.LT.1 )THEN NORM = ZERO ELSE IF( N.EQ.1 )THEN NORM = ABS( X( 1 ) ) ELSE SCALE = ZERO SSQ = ONE * The following loop is equivalent to this call to the LAPACK * auxiliary routine: * CALL DLASSQ( N, X, INCX, SCALE, SSQ ) * DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX IF( X( IX ).NE.ZERO )THEN ABSXI = ABS( X( IX ) ) IF( SCALE.LT.ABSXI )THEN SSQ = ONE + SSQ*( SCALE/ABSXI )**2 SCALE = ABSXI ELSE SSQ = SSQ + ( ABSXI/SCALE )**2 END IF END IF 10 CONTINUE NORM = SCALE * SQRT( SSQ ) END IF * DNRM2 = NORM RETURN * * End of DNRM2. * END ************************************************************************ subroutine dscal(n,da,dx,incx) c c scales a vector by a constant. c uses unrolled loops for increment equal to one. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision da,dx(*) integer i,incx,m,mp1,n,nincx c if( n.le.0 .or. incx.le.0 )return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c nincx = n*incx do 10 i = 1,nincx,incx dx(i) = da*dx(i) 10 continue return c c code for increment equal to 1 c c c clean-up loop c 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end ************************************************************************ integer function idamax(n,dx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c modified 3/93 to return if incx .le. 0. c modified 12/3/93, array(1) declarations changed to array(*) c double precision dx(*),dmax integer i,incx,ix,n c idamax = 0 if( n.lt.1 .or. incx.le.0 ) return idamax = 1 if(n.eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i = 2,n if(dabs(dx(ix)).le.dmax) go to 5 idamax = i dmax = dabs(dx(ix)) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c 20 dmax = dabs(dx(1)) do 30 i = 2,n if(dabs(dx(i)).le.dmax) go to 30 idamax = i dmax = dabs(dx(i)) 30 continue return end ************************************************************************ * Level 2 ************************************************************************ 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*1 TRANS * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' y := alpha*A*x + beta*y. * * TRANS = 'T' or 't' y := alpha*A'*x + beta*y. * * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, m ). * Unchanged on exit. * * X - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - DOUBLE PRECISION array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, the incremented array Y * must contain the vector y. On exit, Y is overwritten by the * updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY * .. External Functions .. LOGICAL LSAME c EXTERNAL LSAME * .. External Subroutines .. c EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * 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 IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF ELSE * * Form y := alpha*A'*x + y. * JY = KY IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = ZERO DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE DO 120, J = 1, N TEMP = ZERO IX = KX DO 110, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * End of DGEMV . * END ************************************************************************ SUBROUTINE DSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, N CHARACTER*1 UPLO * .. Array Arguments .. DOUBLE PRECISION AP( * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DSPMV performs the matrix-vector operation * * y := alpha*A*x + beta*y, * * where alpha and beta are scalars, x and y are n element vectors and * A is an n by n symmetric matrix, supplied in packed form. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the matrix A is supplied in the packed * array AP as follows: * * UPLO = 'U' or 'u' The upper triangular part of A is * supplied in AP. * * UPLO = 'L' or 'l' The lower triangular part of A is * supplied in AP. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * AP - DOUBLE PRECISION array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U' or 'u', the array AP must * contain the upper triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) * and a( 2, 2 ) respectively, and so on. * Before entry with UPLO = 'L' or 'l', the array AP must * contain the lower triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) * and a( 3, 1 ) respectively, and so on. * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. On exit, Y is overwritten by the updated * vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Executable Statements .. * * 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 = 6 ELSE IF( INCY.EQ.0 )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSPMV ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * Set up the start points in X and Y. * IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( N - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( N - 1 )*INCY END IF * * Start the operations. In this version the elements of the array AP * are accessed sequentially with one pass through AP. * * First form y := beta*y. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, N Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, N Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, N Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, N Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO ) $ RETURN KK = 1 IF( LSAME( UPLO, 'U' ) )THEN * * Form y when AP contains the upper triangle. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 60, J = 1, N TEMP1 = ALPHA*X( J ) TEMP2 = ZERO K = KK DO 50, I = 1, J - 1 Y( I ) = Y( I ) + TEMP1*AP( K ) TEMP2 = TEMP2 + AP( K )*X( I ) K = K + 1 50 CONTINUE Y( J ) = Y( J ) + TEMP1*AP( KK + J - 1 ) + ALPHA*TEMP2 KK = KK + J 60 CONTINUE ELSE JX = KX JY = KY DO 80, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO IX = KX IY = KY DO 70, K = KK, KK + J - 2 Y( IY ) = Y( IY ) + TEMP1*AP( K ) TEMP2 = TEMP2 + AP( K )*X( IX ) IX = IX + INCX IY = IY + INCY 70 CONTINUE Y( JY ) = Y( JY ) + TEMP1*AP( KK + J - 1 ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY KK = KK + J 80 CONTINUE END IF ELSE * * Form y when AP contains the lower triangle. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 100, J = 1, N TEMP1 = ALPHA*X( J ) TEMP2 = ZERO Y( J ) = Y( J ) + TEMP1*AP( KK ) K = KK + 1 DO 90, I = J + 1, N Y( I ) = Y( I ) + TEMP1*AP( K ) TEMP2 = TEMP2 + AP( K )*X( I ) K = K + 1 90 CONTINUE Y( J ) = Y( J ) + ALPHA*TEMP2 KK = KK + ( N - J + 1 ) 100 CONTINUE ELSE JX = KX JY = KY DO 120, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO Y( JY ) = Y( JY ) + TEMP1*AP( KK ) IX = JX IY = JY DO 110, K = KK + 1, KK + N - J IX = IX + INCX IY = IY + INCY Y( IY ) = Y( IY ) + TEMP1*AP( K ) TEMP2 = TEMP2 + AP( K )*X( IX ) 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY KK = KK + ( N - J + 1 ) 120 CONTINUE END IF END IF * RETURN * * End of DSPMV . * END ************************************************************************ SUBROUTINE DSPR ( UPLO, N, ALPHA, X, INCX, AP ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX, N CHARACTER*1 UPLO * .. Array Arguments .. DOUBLE PRECISION AP( * ), X( * ) * .. * * Purpose * ======= * * DSPR performs the symmetric rank 1 operation * * A := alpha*x*x' + A, * * where alpha is a real scalar, x is an n element vector and A is an * n by n symmetric matrix, supplied in packed form. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the matrix A is supplied in the packed * array AP as follows: * * UPLO = 'U' or 'u' The upper triangular part of A is * supplied in AP. * * UPLO = 'L' or 'l' The lower triangular part of A is * supplied in AP. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * AP - DOUBLE PRECISION array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U' or 'u', the array AP must * contain the upper triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) * and a( 2, 2 ) respectively, and so on. On exit, the array * AP is overwritten by the upper triangular part of the * updated matrix. * Before entry with UPLO = 'L' or 'l', the array AP must * contain the lower triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) * and a( 3, 1 ) respectively, and so on. On exit, the array * AP is overwritten by the lower triangular part of the * updated matrix. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, K, KK, KX * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Executable Statements .. * * 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 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSPR ', 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 the array AP * are accessed sequentially with one pass through AP. * KK = 1 IF( LSAME( UPLO, 'U' ) )THEN * * Form A when upper triangle is stored in AP. * IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = ALPHA*X( J ) K = KK DO 10, I = 1, J AP( K ) = AP( K ) + X( I )*TEMP K = K + 1 10 CONTINUE END IF KK = KK + J 20 CONTINUE ELSE JX = KX DO 40, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IX = KX DO 30, K = KK, KK + J - 1 AP( K ) = AP( K ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE END IF JX = JX + INCX KK = KK + J 40 CONTINUE END IF ELSE * * Form A when lower triangle is stored in AP. * IF( INCX.EQ.1 )THEN DO 60, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = ALPHA*X( J ) K = KK DO 50, I = J, N AP( K ) = AP( K ) + X( I )*TEMP K = K + 1 50 CONTINUE END IF KK = KK + N - J + 1 60 CONTINUE ELSE JX = KX DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IX = JX DO 70, K = KK, KK + N - J AP( K ) = AP( K ) + X( IX )*TEMP IX = IX + INCX 70 CONTINUE END IF JX = JX + INCX KK = KK + N - J + 1 80 CONTINUE END IF END IF * RETURN * * End of DSPR . * END ************************************************************************ SUBROUTINE DSPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX, INCY, N CHARACTER*1 UPLO * .. Array Arguments .. DOUBLE PRECISION AP( * ), X( * ), Y( * ) * .. * * Purpose * ======= * * DSPR2 performs the symmetric rank 2 operation * * A := alpha*x*y' + alpha*y*x' + A, * * where alpha is a scalar, x and y are n element vectors and A is an * n by n symmetric matrix, supplied in packed form. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the matrix A is supplied in the packed * array AP as follows: * * UPLO = 'U' or 'u' The upper triangular part of A is * supplied in AP. * * UPLO = 'L' or 'l' The lower triangular part of A is * supplied in AP. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * Y - 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. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. INCY must not be zero. * Unchanged on exit. * * AP - DOUBLE PRECISION array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U' or 'u', the array AP must * contain the upper triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) * and a( 2, 2 ) respectively, and so on. On exit, the array * AP is overwritten by the upper triangular part of the * updated matrix. * Before entry with UPLO = 'L' or 'l', the array AP must * contain the lower triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) * and a( 3, 1 ) respectively, and so on. On exit, the array * AP is overwritten by the lower triangular part of the * updated matrix. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. External Subroutines .. EXTERNAL XERBLA * .. * .. Executable Statements .. * * 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( INCY.EQ.0 )THEN INFO = 7 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSPR2 ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * Set up the start points in X and Y if the increments are not both * unity. * IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( N - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( N - 1 )*INCY END IF JX = KX JY = KY END IF * * Start the operations. In this version the elements of the array AP * are accessed sequentially with one pass through AP. * KK = 1 IF( LSAME( UPLO, 'U' ) )THEN * * Form A when upper triangle is stored in AP. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 20, J = 1, N IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN TEMP1 = ALPHA*Y( J ) TEMP2 = ALPHA*X( J ) K = KK DO 10, I = 1, J AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2 K = K + 1 10 CONTINUE END IF KK = KK + J 20 CONTINUE ELSE DO 40, J = 1, N IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN TEMP1 = ALPHA*Y( JY ) TEMP2 = ALPHA*X( JX ) IX = KX IY = KY DO 30, K = KK, KK + J - 1 AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2 IX = IX + INCX IY = IY + INCY 30 CONTINUE END IF JX = JX + INCX JY = JY + INCY KK = KK + J 40 CONTINUE END IF ELSE * * Form A when lower triangle is stored in AP. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 60, J = 1, N IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN TEMP1 = ALPHA*Y( J ) TEMP2 = ALPHA*X( J ) K = KK DO 50, I = J, N AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2 K = K + 1 50 CONTINUE END IF KK = KK + N - J + 1 60 CONTINUE ELSE DO 80, J = 1, N IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN TEMP1 = ALPHA*Y( JY ) TEMP2 = ALPHA*X( JX ) IX = JX IY = JY DO 70, K = KK, KK + N - J AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2 IX = IX + INCX IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX JY = JY + INCY KK = KK + N - J + 1 80 CONTINUE END IF END IF * RETURN * * End of DSPR2 . * END SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ) * .. * * Purpose * ======= * * DTRSV solves one of the systems of equations * * A*x = b, or A'*x = b, * * where b and x are n element vectors and A is an n by n unit, or * non-unit, upper or lower triangular matrix. * * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the equations to be solved as * follows: * * TRANS = 'N' or 'n' A*x = b. * * TRANS = 'T' or 't' A'*x = b. * * TRANS = 'C' or 'c' A'*x = b. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, n ). * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element right-hand side vector b. On exit, X is overwritten * with the solution vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KX LOGICAL NOUNIT * .. External Functions .. LOGICAL LSAME c EXTERNAL LSAME * .. External Subroutines .. c EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRSV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOUNIT = LSAME( DIAG, 'N' ) * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF( LSAME( TRANS, 'N' ) )THEN * * Form x := inv( A )*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 20, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) TEMP = X( J ) DO 10, I = J - 1, 1, -1 X( I ) = X( I ) - TEMP*A( I, J ) 10 CONTINUE END IF 20 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 40, J = N, 1, -1 IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) TEMP = X( JX ) IX = JX DO 30, I = J - 1, 1, -1 IX = IX - INCX X( IX ) = X( IX ) - TEMP*A( I, J ) 30 CONTINUE END IF JX = JX - INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = 1, N IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) TEMP = X( J ) DO 50, I = J + 1, N X( I ) = X( I ) - TEMP*A( I, J ) 50 CONTINUE END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) TEMP = X( JX ) IX = JX DO 70, I = J + 1, N IX = IX + INCX X( IX ) = X( IX ) - TEMP*A( I, J ) 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF END IF ELSE * * Form x := inv( A' )*x. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = X( J ) DO 90, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( I ) 90 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( J ) = TEMP 100 CONTINUE ELSE JX = KX DO 120, J = 1, N TEMP = X( JX ) IX = KX DO 110, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( JX ) = TEMP JX = JX + INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = N, 1, -1 TEMP = X( J ) DO 130, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( I ) 130 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( J ) = TEMP 140 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 160, J = N, 1, -1 TEMP = X( JX ) IX = KX DO 150, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX - INCX 150 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( JX ) = TEMP JX = JX - INCX 160 CONTINUE END IF END IF END IF * RETURN * * End of DTRSV . * END ************************************************************************ * Level 3 ************************************************************************ SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) * .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX, LDA, N CHARACTER*1 UPLO * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ) * .. * * Purpose * ======= * * DSYR performs the symmetric rank 1 operation * * A := alpha*x*x' + A, * * where alpha is a real scalar, x is an n element vector and A is an * n by n symmetric matrix. * * Parameters * ========== * * UPLO - 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. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - DOUBLE PRECISION array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array A must contain the upper * triangular 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. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * max( 1, n ). * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 22-October-1986. * Jack Dongarra, Argonne National Lab. * Jeremy Du Croz, Nag Central Office. * Sven Hammarling, Nag Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KX * .. External Functions .. LOGICAL LSAME c EXTERNAL LSAME * .. External Subroutines .. c EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. * .. Executable Statements .. * * 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 ************************************************************************ SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, $ BETA, C, LDC ) * .. Scalar Arguments .. CHARACTER*1 UPLO, TRANS INTEGER N, K, LDA, LDC DOUBLE PRECISION ALPHA, BETA * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), C( LDC, * ) * .. * * Purpose * ======= * * DSYRK performs one of the symmetric rank k operations * * C := alpha*A*A' + beta*C, * * or * * C := alpha*A'*A + beta*C, * * where alpha and beta are scalars, C is an n by n symmetric matrix * and A is an n by k matrix in the first case and a k by n matrix * in the second case. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the array C is to be referenced as * follows: * * UPLO = 'U' or 'u' Only the upper triangular part of C * is to be referenced. * * UPLO = 'L' or 'l' Only the lower triangular part of C * is to be referenced. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' C := alpha*A*A' + beta*C. * * TRANS = 'T' or 't' C := alpha*A'*A + beta*C. * * TRANS = 'C' or 'c' C := alpha*A'*A + beta*C. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix C. N must be * at least zero. * Unchanged on exit. * * K - INTEGER. * On entry with TRANS = 'N' or 'n', K specifies the number * of columns of the matrix A, and on entry with * TRANS = 'T' or 't' or 'C' or 'c', K specifies the number * of rows of the matrix A. K must be at least zero. * Unchanged on exit. * * ALPHA - DOUBLE PRECISION. * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is * k when TRANS = 'N' or 'n', and is n otherwise. * Before entry with TRANS = 'N' or 'n', the leading n by k * part of the array A must contain the matrix A, otherwise * the leading k by n part of the array A must contain the * matrix A. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. When TRANS = 'N' or 'n' * then LDA must be at least max( 1, n ), otherwise LDA must * be at least max( 1, k ). * Unchanged on exit. * * BETA - DOUBLE PRECISION. * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). * Before entry with UPLO = 'U' or 'u', the leading n by n * upper triangular part of the array C must contain the upper * triangular part of the symmetric matrix and the strictly * lower triangular part of C is not referenced. On exit, the * upper triangular part of the array C 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 C must contain the lower * triangular part of the symmetric matrix and the strictly * upper triangular part of C is not referenced. On exit, the * lower triangular part of the array C is overwritten by the * lower triangular part of the updated matrix. * * LDC - INTEGER. * On entry, LDC specifies the first dimension of C as declared * in the calling (sub) program. LDC must be at least * max( 1, n ). * Unchanged on exit. * * * Level 3 Blas routine. * * -- Written on 8-February-1989. * Jack Dongarra, Argonne National Laboratory. * Iain Duff, AERE Harwell. * Jeremy Du Croz, Numerical Algorithms Group Ltd. * Sven Hammarling, Numerical Algorithms Group Ltd. * * * .. External Functions .. LOGICAL LSAME C EXTERNAL LSAME * .. External Subroutines .. C EXTERNAL XERBLA * .. Intrinsic Functions .. INTRINSIC MAX * .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA DOUBLE PRECISION TEMP * .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Executable Statements .. * * Test the input parameters. * IF( LSAME( TRANS, 'N' ) )THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME( UPLO, 'U' ) * INFO = 0 IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND. $ ( .NOT.LSAME( TRANS, 'T' ) ).AND. $ ( .NOT.LSAME( TRANS, 'C' ) ) )THEN INFO = 2 ELSE IF( N .LT.0 )THEN INFO = 3 ELSE IF( K .LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 7 ELSE IF( LDC.LT.MAX( 1, N ) )THEN INFO = 10 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSYRK ', INFO ) RETURN END IF * * Quick return if possible. * IF( ( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * And when alpha.eq.zero. * IF( ALPHA.EQ.ZERO )THEN IF( UPPER )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, J C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, J C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF ELSE IF( BETA.EQ.ZERO )THEN DO 60, J = 1, N DO 50, I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80, J = 1, N DO 70, I = J, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF * * Start the operations. * IF( LSAME( TRANS, 'N' ) )THEN * * Form C := alpha*A*A' + beta*C. * IF( UPPER )THEN DO 130, J = 1, N IF( BETA.EQ.ZERO )THEN DO 90, I = 1, J C( I, J ) = ZERO 90 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 100, I = 1, J C( I, J ) = BETA*C( I, J ) 100 CONTINUE END IF DO 120, L = 1, K IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 110, I = 1, J C( I, J ) = C( I, J ) + TEMP*A( I, L ) 110 CONTINUE END IF 120 CONTINUE 130 CONTINUE ELSE DO 180, J = 1, N IF( BETA.EQ.ZERO )THEN DO 140, I = J, N C( I, J ) = ZERO 140 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 150, I = J, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE END IF DO 170, L = 1, K IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 160, I = J, N C( I, J ) = C( I, J ) + TEMP*A( I, L ) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE END IF ELSE * * Form C := alpha*A'*A + beta*C. * IF( UPPER )THEN DO 210, J = 1, N DO 200, I = 1, J TEMP = ZERO DO 190, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 190 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 200 CONTINUE 210 CONTINUE ELSE DO 240, J = 1, N DO 230, I = J, N TEMP = ZERO DO 220, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 220 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 230 CONTINUE 240 CONTINUE END IF END IF * RETURN * * End of DSYRK . * END ************************************************************************