C************************************************************************* C * C SPARSE SVD VIA HYBRID BLOCK LANCZOS METHOD ON A'A MATRIX * C * C * C************************************************************************* C * C (C) COPYRIGHT 1991 * C MICHAEL W. BERRY * C ALL RIGHTS RESERVED * C * C PERMISSION TO COPY ALL OR PART OF ANY OF THIS * C SOFTWARE IS ONLY GRANTED UPON APPROVAL FROM THE * C AUTHOR: MICHAEL W. BERRY, DEPT. OF COMPUTER * C SCIENCE, UNIVERSITY OF TENNESSEE, 107 AYRES HALL, * C KNOXVILLE, TN, 37996-1301 (BERRY@CS.UTK.EDU) * C * C * C************************************************************************* C PROGRAM BLS2 C C**************************************************************** C C.......NMAX IS THE DEFAULT MAX NUMBER OF DESIRED TRIPLETS C PARAMETER (NMAX=10, LDV=82, LDU=374) C C**************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) C REAL*8 V(LDV,NMAX),U(LDU,NMAX),S(NMAX),RES(NMAX),TMPV(LDV) REAL*4 ETIME,TMP2(2) C COMMON /COUNT1/ MXVCOUNT COMMON /COUNT2/ MTXVCOUNT INTEGER MXVCOUNT,MTXVCOUNT C C**************************************************************** C C ----- HARWELL-BOEING COLLECTION DATA STRUCTURE C C**************************************************************** C INTEGER NNZERO, NRHS, MATRIX, RESULTS CHARACTER TITLE*72, KEY*8, TYPE*3, PTRFMT*16, INDFMT*16, 1 VALFMT*20, RHSFMT*20, NAME*30 LOGICAL VECTORS C C----------------------------------------------------------- C USER-SUPPLIED PARAMETERS FOR BLKLAN: C C SET NRMAX = ROW DIMENSION OF MATRIX A C SET NCMAX = (COLUMN DIMENSION OF MATRIX A)+1 C (ASSUMING NO. OF ROWS > NO. OF COLS) C NZMAX = MAXIMUM NO. OF NONZEROS IN MATRIX A C PARAMETER(NRMAX=500,NCMAX=200,NZMAX=2000) REAL*8 VALUE(NZMAX), ZTEMP(NRMAX,NCMAX) INTEGER POINTR(NCMAX), ROWIND(NZMAX), NROW, NCOL COMMON /MATRIXA/ NROW,NCOL,VALUE,POINTR,ROWIND,ZTEMP C----------------------------------------------------------- C MATRIX=7 IN=1 RESULTS=2 IVEC=3 MXVCOUNT=0 MTXVCOUNT=0 C C**************************************************************** C C OPEN FILES FOR INPUT/OUTPUT C OPEN(IN, FILE='BLI1') OPEN(MATRIX, FILE='LAI7') OPEN(RESULTS,FILE='BLO2') C READ(IN,*) NAME, NC, NB, NUMS, TOL, VECTORS IF(VECTORS) OPEN(IVEC,FILE='BLO3',FORM='UNFORMATTED') C C.......GET MATRIX A C READ (MATRIX,10) TITLE, KEY, 1 TYPE, NROW, NCOL, NNZERO, NRHS, PTRFMT, INDFMT, VALFMT, RHSFMT 10 FORMAT (A72, A8 // A3, 11X, 4I14 / 2A16, 2A20) C C LEAVE IF MATRIX IS TOO BIG ---- C IF (NCOL.GT. NCMAX .OR. NNZERO .GT. NZMAX .OR. 1 NROW.GT. NRMAX) THEN WRITE(*,*) ' SORRY, YOUR MATRIX IS TOO BIG! ' STOP ENDIF C C READ DATA... C READ (MATRIX,PTRFMT) (POINTR (I), I = 1, NCOL+1) READ (MATRIX,INDFMT) (ROWIND (I), I = 1, NNZERO) READ (MATRIX,VALFMT) (VALUE (I), I = 1, NNZERO) C C DEFINE LAST ELEMENT OF POINTR IN CASE IT IS NOT. C POINTR(NCOL+1) = NNZERO + 1 C C---------------------------------------------------------------------------- C M=NROW N=NCOL NCOLD=NC NBOLD=NB NK=NUMS NKOLD=NK C C-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= T1=ETIME(TMP2(1)) C CALL BLKLAN(N,NKOLD,LDV,V,S,NCOLD,NBOLD,ITER, * TOL,RES,150,NK,NC,NB,IMEM) C TIME=ETIME(TMP2(1))-T1 C IF(VECTORS) THEN C T1=ETIME(TMP2(1)) DO 35 JJ = 1,NK C C...........MULTIPLY BY MATRIX A'A FIRST C CALL OPN(N,V(1,JJ),TMPV) TMP =DDOT(N,V(1,JJ),1,TMPV,1) CALL DAXPY(N,-TMP,V(1,JJ),1,TMPV,1) TMP =DSQRT(TMP) XNORM=DSQRT(DDOT(N,TMPV,1,TMPV,1)) C C...........MULTIPLY BY MATRIX A TO GET (SCALED) LEFT S-VECTOR C CALL OPA(V(1,JJ),U(1,JJ)) TMP1=1.0D0/TMP CALL DSCAL(M,TMP1,U(1,JJ),1) XNORM=XNORM*TMP1 RES(JJ)=XNORM S(JJ)=TMP C 35 CONTINUE C TIME=TIME+ETIME(TMP2(1))-T1 C ENDIF C C.......WRITE UNFORMATTED FILE CONTAINING LEFT AND RIGHT S-VECTOR PAIRS C IF(VECTORS) THEN DO 510 J = 1,NK WRITE(IVEC) (U(I,J), I=1,M) WRITE(IVEC) (V(I,J), I=1,N) 510 CONTINUE ENDIF C WRITE(RESULTS,2000) ITER,NKOLD,NK,NBOLD,NB,NCOLD,NC, * MXVCOUNT,MTXVCOUNT,MXVCOUNT+MTXVCOUNT,IMEM, * VECTORS,TOL,TITLE,NAME,NROW,NCOL C 2000 FORMAT( * 1X,'... ' * /1X,'... HYBRID BLOCK LANCZOS [A^TA]', * /1X,'... NO. OF ITERATIONS TAKEN =',I10 * /1X,'... NO. OF TRIPLETS SOUGHT =',I10 * /1X,'... NO. OF TRIPLETS FOUND =',I10 * /1X,'... INITIAL BLOCKSIZE =',I10 * /1X,'... FINAL BLOCKSIZE =',I10 * /1X,'... MAXIMUM SUBSPACE BOUND =',I10 * /1X,'... FINAL SUBSPACE BOUND =',I10 * /1X,'... NO. MULTIPLICATIONS BY A =',I10 * /1X,'... NO. MULT. BY TRANSPOSE(A)=',I10 * /1X,'... TOTAL SPARSE MAT-VEC MULT.=',I10 * /1X,'... MEMORY NEEDED IN BYTES =',I10 * /1X,'... WANT S-VECTORS ON OUTPUT? ',L10 * /1X,'... [T/F] ' * /1X,'... TOLERANCE =',1PE10.2 * //1X,A50/1X,A50 * /1X,'... NO. OF DOCS (ROWS OF A) = ',I8 * /1X,'... NO. OF TERMS (COLS OF A) = ',I8/1X,'... ') WRITE(RESULTS,9994) TIME WRITE(RESULTS,9997) (I,S(I),RES(I),I = 1,NK) 9994 FORMAT(/1X,'...... BLSVD EXECUTION TIME=',1PE10.2) 9997 FORMAT(1X,'...... ' * /1X,'...... ',3X,3X,' COMPUTED S-VALUES ',2X, * '(','RES. NORMS',')' * /1X,'...... ' * /(1X,'...... ',I3,3X,1PE22.14,2X,'(',1PE11.2,')')) STOP END SUBROUTINE BLKLAN(N,IK,LDV,V,EIG,IC,IB,ITER,TOL,RES,MAXIT, * IKO,ICO,IBO,IMEM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,LDV,IB,IC,IK,ITER,MAXIT REAL*8 V(LDV,IK),EIG(IK),RES(IK),RANDOM C C *************************************** C * * C * [BLOCK LANCZOS SVD ALGORITHM] * C * * C * DETERMINE IK-LARGEST E-PAIRS OF * C * A'A WHERE A IS LARGE AND SPARSE. * C * * C * FIVE PHASES: * C * * C * PHASE 1: BLOCK LANCZOS TO YIELD * C * SYM. BLK-TRIDIAG. MATRIX S WHICH * C * SHARES THE SAME E-VALUES AS A'A. * C * TOTAL RE-ORTHOGONALIZATION USED * C * * C * PHASE 2: LANCZOS METHOD FOR TRI- * C * DIAGONALIZING THE S MATRIX FROM * C * PHASE 1 TO YIELD THE TRI-DIAG. * C * MATRIX T WHICH PRESERVES THE * C * SAME E-VALUES. * C * * C * PHASE 2A: POINT LANCZOS WITH * C * RE-ORTHOGONALIZATION FOR NB=1 * C * OR IB=1 CASE. * C * * C * PHASE 3: APPLY THE STANDARD QL * C * ITERATION TO DIAGONALIZE T AND * C * HENCE PRODUCE THE APPROX. E-VALUES * C * (ARRAY ALPHA) OF MATRIX A. * C * * C * PHASE 4: CONVERGENCE TEST. * C * * C * PHASE 5: ITERATION RESTART. * C * PROJECT AGAINST CONVERGED * C * E-VECTORS BEFORE PHASE 1 STARTS. * C * * C *************************************** C * * C * INPUT: IC = MAX(IB*IS), * C * WHERE IB = BLOCKSIZE, * C * IS = NUMBER OF BLOCKS * C * IK = NO. DESIRED E-PAIRS * C * TOL = RESIDUAL TOLERANCE * C * MAXIT = MAXIMUM # ITERATIONS * C * * C * OUTPUT: IT = NO. OF ITERATIONS * C * RES(IKO) = RESIDUALS * C * IBO = FINAL BLOCKSIZE * C * ICO = FINAL MEM BOUND * C * IMEM = STORAGE USED * C * (IN BYTES) * C * IKO = NO. E-PAIRS FOUND * C * EIG(IKO) = E-VALUES * C * (DESCENDING ORDER) * C * V(LDV,IKO) = CORRESP. E-VECTORS * C * * C * REQUIRED: * C * * C * CALLS OPN(N,X,Y) FOR Y=A'A*X * C * CALLS OPM(N,NC,X,LDX,Y,LDY) * C * FOR Y=A'A*X (BLOCK OF VECTORS) * C * * C * EXTERNALS ROUTINES: * C * * C * RANDOM - RANDOM NO. GEN. * C * DDOT - DOTPRODUCT (BLAS1) * C * DNRM2 - E-NORM (BLAS1) * C * DAXPY - Y=Y+A*X (BLAS1) * C * DGEMM - C=C-A*B (LAPACK) * C * DGEMV - Y=Y+A*X (LAPACK) * C * ORTHG - MGS WITH PROJECTION * C * DSBMV - BAND MATRIX-VECTOR MULT.* C * (LAPACK) * C * TQL2 - EISPACK QL METHOD * C *************************************** C C----------------------------------------------------------- C USER-SUPPLIED PARAMETERS FOR BLKLAN: C C SET NRMAX = ROW DIMENSION OF MATRIX A C SET NCMAX = (COLUMN DIMENSION OF MATRIX A)+1 C (ASSUMING NO. OF ROWS > NO. OF COLS) C NZMAX = MAXIMUM NO. OF NONZEROS IN MATRIX A C PARAMETER(NRMAX=500,NCMAX=200,NZMAX=2000) REAL*8 VALUE(NZMAX), ZTEMP(NRMAX,NCMAX) INTEGER POINTR(NCMAX), ROWIND(NZMAX), NROW, NCOL COMMON /MATRIXA/ NROW,NCOL,VALUE,POINTR,ROWIND,ZTEMP C----------------------------------------------------------- C COMMON /COUNT1/ MXVCOUNT COMMON /COUNT2/ MTXVCOUNT INTEGER MXVCOUNT,MTXVCOUNT C C AUXILIARY STORAGE DECLARATIONS: C C--------------------------------------------------------- C USER-SUPPLIED PARAMETERS FOR BLKLAN: C C SET IB2 = INITIAL BLOCK SIZE (NB) C IK2 = NO. DESIRED S-VALUES (NUMS) C IC2 = INITIAL SUBSPACE DIMENSION (NC) C N2 = COL DIMENSION OF MATRIX A (N) C--------------------------------------------------------- C PARAMETER( * IB2=10, * IK2=10, * IC2=50, * N2=82) C REAL*8 W(N2,IB2),Y(N2,IB2),S(IB2,IC2),R(IB2,IC2-IB2), * TEMP(N2,IC2), VV(N2,IC2),V0(N2,IK2), TRES(IB2), * UVTMP(N2,IC2+IK2), BIGS(IB2+1,IC2),P(N2),Q(N2),T(N2), * PP(IC2,IC2),QQ(IC2,IC2),ALPHA(IC2), Z(N2),BETA(IC2) C INTEGER KINDX(IB2) C C--------------------------------------------------------- C C C NS IS THE NUMBER OF BLOCK COLUMN PARTITIONS TO FORM C C NB IS THE NUMBER OF COLUMNS IN V(J) C C I.E., V = [ V(1), V(2), ... , V(NS) ] C C......SET DIMENSIONAL PARAMETERS (LOCAL) C NB=IB NC=IC K=IK C C----------------------------------------------------------------------- C IMEM=LDV*IK2+IK2+IK2+ * N2*IB2+N2*IB2+IB2*IC2+IB2*(IC2-IB2)+ * N2*IK2+N2*IC2+N2*IC2+ * N2*(IC2+IK2)+IK2+ * IK2+(IB2+1)*IC2+N2+2*N2+ * N2+IC2+2*IC2*IC2+IC2+IB2 C IMEM=8*IMEM C C......INTEGER STORAGE C IMEM=IMEM+4*IB2 C C......IMEM STORES TOTAL AUXILIARY STORAGE REQUIRED (IN BYTES) C C----------------------------------------------------------------------- C C......ITERATION COUNTER: ITER C C......CONVERGED VECTOR COUNTER : ICONV C LDB=IB2+1 LDR=IB2 LDS=IB2 ITER=0 ICONV=0 K0=0 NBOLD=0 C C ******************** B E G I N P H A S E 1 ********************** C C......DETERMINE NUMBER OF BLOCKS FOR 1ST ITERATION (NS) C NS=NC/NB IF(NS.LT.2) THEN NB=NC/2 NS=NC/NB ENDIF C 1000 ITER=ITER+1 C NN=NB*NS C IF(NB.EQ.1) THEN DO 195 II=1,N P(II)=VV(II,1) 195 CONTINUE C C.........RESET CONVERGED VECTORS COUNTER C IF(NBOLD.NE.1) K0=0 GO TO 888 C ENDIF C IF(ITER.EQ.1) THEN C C.........CHOOSE AN INITIAL (RANDOM) V[1] MATRIX (N BY NB) C IRAND=91827211 DO 10 J=1,NB DO 20 I=1,N VV(I,J)=RANDOM(IRAND) 20 CONTINUE 10 CONTINUE C CALL ORTHG(N,0,NB,Y,N2,VV,N2,TEMP) C ENDIF C C C......INITIALIZE S AND R ARRAYS C C DO 315 JJ=1,NB*NS DO 316 II=1,NB S(II,JJ)=0.0D0 316 CONTINUE 315 CONTINUE C IF(NS.GT.1) THEN DO 415 JJ=1,NB*NS-NB DO 416 II=1,NB R(II,JJ)=0.0D0 416 CONTINUE 415 CONTINUE ENDIF C C COMPUTE S(1) = V(1)'*A*V(1) C CALL OPM(N,NB,VV(1,1),N2,Y(1,1),N2) C C LOAD S(1) SUBMATRIX C CALL DGEMM('T','N',NB,NB,N,1.0D0,VV,N2,Y,N2,0.0D0,S(1,1),LDS) C C ITERATE FOR J = 2, ... , NS C DO 30 J =2,NS C C COMPUTE POINTER TO CURRENT BLOCK PARTITION C JPTR=(J-2)*NB+1 C C C----------------------------------------------------------------------- C COMPUTE Y(J) = A * V(J-1) - V(J-1) * S(J-1) - V(J-2) * R(J-2)' C----------------------------------------------------------------------- C CALL DGEMM('N','N',N,NB,NB,-1.0D0,VV(1,JPTR),N2, * S(1,JPTR),LDS,1.0D0,Y,N2) C IF(J.GT.2) THEN C CALL DGEMM('N','T',N,NB,NB,-1.0D0,VV(1,JPTR-NB),N2, * R(1,JPTR-NB),LDR,1.0D0,Y,N2) ENDIF C C..........SAVE RESIDUALS FOR CONVERGENCE TEST IN PHASE 4 C (STORED IN Y(1)) C IF(J.EQ.2.AND.ITER.GT.1) THEN C DO 65 JJ=1,NB TRES(JJ)=DNRM2(N,Y(1,JJ),1) 65 CONTINUE C ENDIF C C-------------------------------------------------------------- C..........ORTHOGONALIZE Y(J) W.R.T. [V(1),....V(J) | V(0) ] C-------------------------------------------------------------- C IF(ICONV.GT.0) THEN C C...............ORTHOGONALIZE Y(J) W.R.T. [V(0), V(1), ... , V(J)] C NJ=NB*(J-1)+ICONV DO 731 KK=1,NJ DO 732 LL=1,N IF(KK.LE.NB*(J-1)) THEN UVTMP(LL,KK)=VV(LL,KK) ELSE UVTMP(LL,KK)=V0(LL,KK-NB*(J-1)) ENDIF 732 CONTINUE 731 CONTINUE C ELSE C C...............ORTHOGONALIZE Y(J) W.R.T. [V(1), ... ,V(J)] C NJ=NB*(J-1) DO 131 KK=1,NJ DO 132 LL=1,N UVTMP(LL,KK)=VV(LL,KK) 132 CONTINUE 131 CONTINUE ENDIF DO 431 KK=1,NB DO 432 LL=1,N UVTMP(LL,NJ+KK)=Y(LL,KK) Y(LL,KK)=0.0D0 432 CONTINUE 431 CONTINUE C C---------------------------------------------- CALL ORTHG(N,NJ,NB,Y,N2,UVTMP,N2,TEMP) C---------------------------------------------- C C GET QR FACTORIZATION: Y(J) = V(J) * R(J-1), [ R(J-1) := UPPER TRI.] C C INCREMENT POINTER FOR V(J) C JINC=JPTR+NB C C V(J) SUBMATRIX NOW LOADED C C LOAD R(J) SUBMATRIX C DO 100 JJ = 1,NB DO 110 I= 1,JJ R(I,JJ+JPTR-1)=Y(I,JJ) 110 CONTINUE 100 CONTINUE C C LOAD S(J+1) SUBMATRIX C DO 105 JJ = 1,NB DO 115 II= 1,N VV(II,JINC+JJ-1)=UVTMP(II,NJ+JJ) 115 CONTINUE 105 CONTINUE C CALL OPM(N,NB,VV(1,JINC),N2,Y(1,1),N2) C CALL DGEMM('T','N',NB,NB,N,1.0D0,VV(1,JINC),N2, * Y,N2,0.0D0,S(1,JINC),LDS) C C 30 CONTINUE C C----------------------------------------------------------------- C FORM THE SYMMETRIC BLOCK TRI-DIAGONAL MATRIX S IN BIGS ARRAY C----------------------------------------------------------------- C CALL FORMBIGS(NN,NB,R,S,BIGS,LDB,LDR,LDS) C C ******************** B E G I N P H A S E 2 ********************** C C LANCZOS TRI-DIAGONALIZATION: PP GENERATES E-VECTORS OF BAND MATRIX C C CHOOSE A STARTING VECTOR P(1) [HAVING UNITY 2-NORM] C IHBW=NB C DO 200 I=1,NN P(I)=RANDOM(IRAND) 200 CONTINUE PNM=DNRM2(NN,P,1) DO 201 I=1,NN P(I)=P(I)/PNM 201 CONTINUE C DO 206 II = 1,NN PP(II,1)=P(II) Z(II)=0.0D0 206 CONTINUE C C ITERATE FOR J = 1, 2, ... , NN C C DO 210 J = 1,NN C CALL DSBMV('U',NN,IHBW,1.0D0,BIGS,LDB,P,1,1.0D0,Z,1) ALPHA(J)=DDOT(NN,P,1,Z,1) IF(J.EQ.NN) GO TO 500 C C COMPUTE Z(J) := Z(J)-ALPHA(J)*P C CALL DAXPY(NN,-ALPHA(J),P,1,Z,1) C C----------------------------------------------------------------------- C.......ORTHOGONALIZE Z W.R.T. PREVIOUS PP'S C----------------------------------------------------------------------- C DO 1620 JJ=1,J DO 1625 II=1,NN UVTMP(II,JJ)=PP(II,JJ) 1625 CONTINUE 1620 CONTINUE IF(J.GT.1) THEN DO 1626 II=1,NN UVTMP(II,J+1)=Z(II) 1626 CONTINUE CALL ORTHG(NN,J,1,W,N2,UVTMP,N2,TEMP) DO 1627 II=1,NN Z(II)=UVTMP(II,J+1) 1627 CONTINUE ELSE DUM=-DDOT(NN,UVTMP(1,1),1,Z,1) CALL DAXPY(NN,DUM,UVTMP(1,1),1,Z,1) ENDIF C C COMPUTE BETA(J) C BETA(J)=DNRM2(NN,Z,1) C IF(BETA(J).NE.0.0D0) THEN C C COMPUTE P(J+1) := Z(J)/BETA(J) C DO 230 JJ = 1,NN T(JJ)=P(JJ) 230 CONTINUE DO 231 JJ = 1,NN P(JJ)=Z(JJ)/BETA(J) 231 CONTINUE DO 235 II = 1,NN PP(II,J+1)=P(II) Z(II)=-BETA(J)*T(II) 235 CONTINUE ELSE GO TO 500 ENDIF C 210 CONTINUE C C C ******************** B E G I N P H A S E 2A (NB=1 CASE) ********** C C (POINT ALGORITHM) C C LANCZOS TRI-DIAGONALIZATION: VV GENERATES E-VECTORS (VL=N) C C C C----------------------------------------------------------------- C C ITERATE FOR J =1, ... , NN C 888 CONTINUE DO 309 JJ=1,N Z(JJ)=0.0D0 309 CONTINUE DO 310 J = 1,NN C CALL OPN(N,P,Q) CALL DAXPY(N,1.0D0,Q,1,Z,1) ALPHA(J)=DDOT(N,P,1,Z,1) IF(ALPHA(J).NE.0.0D0) THEN C C IF(J.EQ.NN) GO TO 500 C C COMPUTE Z(J) := Z(J)-ALPHA(J)*P C CALL DAXPY(N,-ALPHA(J),P,1,Z,1) C IF(J.EQ.1) THEN ZNM=DNRM2(N,Z,1) IF(ZNM.LE.TOL) THEN DO 327 II=1,N V(II,ICONV+1)=P(II) 327 CONTINUE EIG(ICONV+1)=ALPHA(1) RES(ICONV+1)=ZNM K=K-1 ITER=ITER-1 GO TO 999 ENDIF ENDIF C C----------------------------------------------------------------------- C.......ORTHOGONALIZE Z W.R.T. CONVERGED RIGHT S-VECTORS AND PREVIOUS VV'S C----------------------------------------------------------------------- C DO 620 JJ=1,ICONV+J DO 625 II=1,N IF(JJ.LE.ICONV) THEN UVTMP(II,JJ)=V0(II,JJ) ELSE UVTMP(II,JJ)=VV(II,JJ-ICONV) ENDIF 625 CONTINUE 620 CONTINUE C IF(ICONV+J.GT.1) THEN DO 1921 II=1,N UVTMP(II,ICONV+J+1)=Z(II) 1921 CONTINUE CALL ORTHG(N,ICONV+J,1,W,N2,UVTMP,N2,TEMP) DO 1922 II=1,N Z(II)=UVTMP(II,ICONV+J+1) 1922 CONTINUE ELSE DUM=-DDOT(N,UVTMP(1,1),1,Z,1) CALL DAXPY(N,DUM,UVTMP(1,1),1,Z,1) ENDIF C C----------------------------------------------------------------- C C COMPUTE BETA(J) C BETA(J)=DNRM2(N,Z,1) C IF(BETA(J).NE.0.0D0) THEN C C COMPUTE P(J+1) := Z(J)/BETA(J) C DO 330 JJ = 1,N T(JJ)=P(JJ) 330 CONTINUE DO 331 JJ = 1,N P(JJ)=Z(JJ)/BETA(J) 331 CONTINUE DO 332 JJ = 1,N VV(JJ,J+1)=P(JJ) Z(JJ)=-BETA(J)*T(JJ) 332 CONTINUE C C----------------------------------------------------------------- C ELSE GO TO 500 ENDIF ELSE GO TO 500 ENDIF C 310 CONTINUE C 500 CONTINUE C C ******************** B E G I N P H A S E 3 ********************** C C SOLVE SYMMETRIC TRIDIAGONAL EVP C DO 735 JJ=1,NN DO 736 II=1,NN QQ(II,JJ)=0.0D0 736 CONTINUE 735 CONTINUE DO 737 II=1,NN QQ(II,II)=1.0D0 737 CONTINUE C DO 738 JJ=NN,2,-1 BETA(JJ)=BETA(JJ-1) 738 CONTINUE BETA(1)=0.0D0 C CALL TQL2(IC2,NN,ALPHA,BETA,QQ,INFO) C C......RESORT ALPHA'S AND COLUMNS OF PP (NEED DESCENDING ORDER) C DO 740 JJ = 1,NN Z(JJ)=ALPHA(NN-JJ+1) DO 741 II=1,NN UVTMP(II,JJ)=QQ(II,NN-JJ+1) 741 CONTINUE 740 CONTINUE C DO 742 JJ = 1,NN ALPHA(JJ)=Z(JJ) DO 743 II=1,NN QQ(II,JJ)=UVTMP(II,JJ) 743 CONTINUE 742 CONTINUE C C ******************** B E G I N P H A S E 4 ********************** C NBOLD=NB IF(ITER.GT.1.AND.NB.GT.1) THEN C K0=0 DO 410 JJ=1,MIN0(NB,K) IF(DABS(TRES(JJ)).LE.TOL) THEN K0=K0+1 KINDX(K0)=JJ EIG(ICONV+K0)=ALPHA(JJ) RES(ICONV+K0)=TRES(JJ) ENDIF 410 CONTINUE C IF(NB.GE.K) THEN NB=NB-K0 ELSE NB=MIN0(NB,K-K0) ENDIF C NC=NC-K0 K=K-K0 IF(K.LE.0) GO TO 997 C C......DETERMINE NUMBER OF BLOCKS FOR NEXT ITERATION (NS) C NS=NC/NB IF(NS.LT.2) THEN NB=NC/2 NS=NC/NB ENDIF ENDIF C C ******************** B E G I N P H A S E 5 ********************** C 997 IF(NBOLD.GT.1) THEN CALL DGEMM('N','N',NN,NN,NN,1.0D0,PP,IC2,QQ,IC2, * 0.0D0,TEMP,N2) CALL DGEMM('N','N',N,NB+K0,NN,1.0D0,VV,N2, * TEMP,N2,0.0D0,V,LDV) ELSE CALL DGEMV('N',N,NN,1.0D0,VV,N2,QQ(1,1),1,0.0D0,V(1,1),1) ENDIF C IF(K0.GT.0) THEN DO 420 JJ=ICONV+1,ICONV+K0 DO 435 II=1,N V0(II,JJ)=V(II,KINDX(JJ-ICONV)) 435 CONTINUE 420 CONTINUE C ICONV=ICONV+K0 IF(K.LE.0) THEN ITER=ITER-1 GO TO 999 ENDIF C ENDIF C C......RELOAD UNCONVERGED E-VECTORS C IF(K0.GT.0) THEN C KK=0 DO 440 JJ=1,NB+K0 DO 444 LL=1,K0 IF(KINDX(LL).EQ.JJ) GO TO 440 444 CONTINUE KK=KK+1 DO 445 II=1,N VV(II,KK)=V(II,JJ) 445 CONTINUE 440 CONTINUE ELSE DO 640 JJ=1,NB DO 645 II=1,N VV(II,JJ)=V(II,JJ) 645 CONTINUE 640 CONTINUE C ENDIF C IF(ITER.EQ.MAXIT) GO TO 999 C GO TO 1000 C C......RETURN ALL CONVERGED S-VECTORS (RIGHT AND LEFT) C 999 DO 920 JJ=1,ICONV DO 935 II=1,N V(II,JJ)=V0(II,JJ) 935 CONTINUE 920 CONTINUE C IKO=IK-K IBO=NB ICO=NC C RETURN END C SUBROUTINE FORMBIGS(N,NDP,R,S,BIGS,LDB,LDR,LDS) IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,LDB,LDR,LDS REAL*8 S(LDS,N),R(LDR,N-NDP),BIGS(LDB,N) C C ************************************** C * * C * ROUTINE TO FORM THE BLOCK UPPER- * C * TRIANGULAR MATRIX S FROM THE * C * BLOCK LANCZOS ALGORITHM IN PHASE * C * 1 OF ROUTINE BLKLAN. ARRAY BIGS * C * WILL CONTAIN THIS MATRIX IN THE * C * FOLLOWING BAND MATRIX FORMAT: * C * * C * EXAMPLE WITH 4 SUP. DIAGONALS: * C * * C * [0 0 0 0--1ST SUP. DIAGONAL----] * C * [0 0 0 ---2ND SUP. DIAGONAL----] * C * [0 0 -----3RD SUP. DIAGONAL----] * C * [0 -------4TH SUP. DIAGONAL----] * C * [-------- MAIN DIAGONAL--------] * C * * C * THE FOLLOWING DATA STRUCTURE * C * IS ASSUMED FOR THE UPPER-TRI. * C * SUBMATRICES S(J) AND R(J), WHERE * C * J = 1,2,...P. [HERE P IS THE * C * NO. OF BLOCK PARTITIONS USED IN * C * BLKLAN]. * C * * C * S := [S(1) | S(2) | ... S(P)] * C * R := [R(1) | R(2) | ... R(P-1)] * C * * C * BIGS IS N BY N. * C * S(J) IS NDP BY NDP. * C * R(J) IS NDP BY NDP. * C * * C * NOTE: NDP := N/P * C * * C ************************************** C C NN=NDP+1 C C LOAD MAIN DIAGONAL OF BIGS C DO 10 I = 1,N,NDP DO 20 K = 1,NDP BIGS(NN,I+K-1)=S(K,I+K-1) 20 CONTINUE 10 CONTINUE C C LOAD SUPER DIAGONALS OF BIGS (TOP TO BOTTOM ORDER) C DO 30 I = 1,NDP C C PAD ZEROS ON THE LEFT C DO 40 K = 1,NDP-I+1 BIGS(I,K)=0.0D0 40 CONTINUE C IF(I.EQ.1) THEN C IF(NDP.NE.N) THEN DO 50 K=NDP-I+2,N,NDP KINC=((K/NDP)-1)*NDP DO 60 KK=1,NDP BIGS(I,K+KK-1)=R(KK,KK+KINC) 60 CONTINUE 50 CONTINUE ENDIF C ELSE C DO 70 K=NDP-I+2,N,NDP IF(I.NE.2) THEN KINC=(K/NDP)*NDP ELSE KINC=((K/NDP)-1)*NDP ENDIF C C LOAD ELEMENTS FROM S(J) SUBMATRICES C DO 80 KK=1,I-1 BIGS(I,K+KK-1)=S(KK,NDP-I+KK+1+KINC) 80 CONTINUE C C LOAD ELEMENTS FROM R^T(J) SUBMATRICES C DO 90 KK=1,NDP-I+1 BIGS(I,K+I-1+KK-1)=R(KK,I+KK-1+KINC) 90 CONTINUE C 70 CONTINUE C ENDIF C 30 CONTINUE C RETURN END SUBROUTINE DSBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. SCALAR ARGUMENTS .. REAL*8 ALPHA, BETA INTEGER INCX, INCY, K, LDA, N CHARACTER*1 UPLO * .. ARRAY ARGUMENTS .. REAL*8 A( LDA, * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * DSBMV 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 BAND MATRIX, WITH K SUPER-DIAGONALS. * * PARAMETERS * ========== * * UPLO - CHARACTER*1. * ON ENTRY, UPLO SPECIFIES WHETHER THE UPPER OR LOWER * TRIANGULAR PART OF THE BAND MATRIX A IS BEING SUPPLIED AS * FOLLOWS: * * UPLO = 'U' OR 'U' THE UPPER TRIANGULAR PART OF A IS * BEING SUPPLIED. * * UPLO = 'L' OR 'L' THE LOWER TRIANGULAR PART OF A IS * BEING SUPPLIED. * * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE ORDER OF THE MATRIX A. * N MUST BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * K - INTEGER. * ON ENTRY, K SPECIFIES THE NUMBER OF SUPER-DIAGONALS OF THE * MATRIX A. K MUST SATISFY 0 .LE. K. * UNCHANGED ON EXIT. * * ALPHA - REAL*8. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * A - REAL*8 ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY WITH UPLO = 'U' OR 'U', THE LEADING ( K + 1 ) * BY N PART OF THE ARRAY A MUST CONTAIN THE UPPER TRIANGULAR * BAND PART OF THE SYMMETRIC MATRIX, SUPPLIED COLUMN BY * COLUMN, WITH THE LEADING DIAGONAL OF THE MATRIX IN ROW * ( K + 1 ) OF THE ARRAY, THE FIRST SUPER-DIAGONAL STARTING AT * POSITION 2 IN ROW K, AND SO ON. THE TOP LEFT K BY K TRIANGLE * OF THE ARRAY A IS NOT REFERENCED. * THE FOLLOWING PROGRAM SEGMENT WILL TRANSFER THE UPPER * TRIANGULAR PART OF A SYMMETRIC BAND MATRIX FROM CONVENTIONAL * FULL MATRIX STORAGE TO BAND STORAGE: * * DO 20, J = 1, N * M = K + 1 - J * DO 10, I = MAX( 1, J - K ), J * A( M + I, J ) = MATRIX( I, J ) * 10 CONTINUE * 20 CONTINUE * * BEFORE ENTRY WITH UPLO = 'L' OR 'L', THE LEADING ( K + 1 ) * BY N PART OF THE ARRAY A MUST CONTAIN THE LOWER TRIANGULAR * BAND PART OF THE SYMMETRIC MATRIX, SUPPLIED COLUMN BY * COLUMN, WITH THE LEADING DIAGONAL OF THE MATRIX IN ROW 1 OF * THE ARRAY, THE FIRST SUB-DIAGONAL STARTING AT POSITION 1 IN * ROW 2, AND SO ON. THE BOTTOM RIGHT K BY K TRIANGLE OF THE * ARRAY A IS NOT REFERENCED. * THE FOLLOWING PROGRAM SEGMENT WILL TRANSFER THE LOWER * TRIANGULAR PART OF A SYMMETRIC BAND MATRIX FROM CONVENTIONAL * FULL MATRIX STORAGE TO BAND STORAGE: * * DO 20, J = 1, N * M = 1 - J * DO 10, I = J, MIN( N, J + K ) * A( M + I, J ) = MATRIX( I, J ) * 10 CONTINUE * 20 CONTINUE * * 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 * ( K + 1 ). * UNCHANGED ON EXIT. * * X - REAL*8 ARRAY OF DIMENSION AT LEAST * ( 1 + ( N - 1 )*ABS( INCX ) ). * 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 - REAL*8. * ON ENTRY, BETA SPECIFIES THE SCALAR BETA. * UNCHANGED ON EXIT. * * Y - REAL*8 ARRAY OF DIMENSION AT LEAST * ( 1 + ( N - 1 )*ABS( INCY ) ). * BEFORE ENTRY, THE INCREMENTED ARRAY Y MUST CONTAIN THE * 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 .. REAL*8 ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. REAL*8 TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, KPLUS1, KX, KY, L * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX, MIN * .. * .. 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( K.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.( K + 1 ) )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( 'DSBMV ', 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 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, 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 IF( LSAME( UPLO, 'U' ) )THEN * * FORM Y WHEN UPPER TRIANGLE OF A IS STORED. * KPLUS1 = K + 1 IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 60, J = 1, N TEMP1 = ALPHA*X( J ) TEMP2 = ZERO L = KPLUS1 - J DO 50, I = MAX( 1, J - K ), J - 1 Y( I ) = Y( I ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + A( L + I, J )*X( I ) 50 CONTINUE Y( J ) = Y( J ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2 60 CONTINUE ELSE JX = KX JY = KY DO 80, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO IX = KX IY = KY L = KPLUS1 - J DO 70, I = MAX( 1, J - K ), J - 1 Y( IY ) = Y( IY ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + A( L + I, J )*X( IX ) IX = IX + INCX IY = IY + INCY 70 CONTINUE Y( JY ) = Y( JY ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY IF( J.GT.K )THEN KX = KX + INCX KY = KY + INCY END IF 80 CONTINUE END IF ELSE * * FORM Y WHEN LOWER TRIANGLE OF A IS STORED. * 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*A( 1, J ) L = 1 - J DO 90, I = J + 1, MIN( N, J + K ) Y( I ) = Y( I ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + A( L + I, J )*X( I ) 90 CONTINUE Y( J ) = Y( J ) + ALPHA*TEMP2 100 CONTINUE ELSE JX = KX JY = KY DO 120, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO Y( JY ) = Y( JY ) + TEMP1*A( 1, J ) L = 1 - J IX = JX IY = JY DO 110, I = J + 1, MIN( N, J + K ) IX = IX + INCX IY = IY + INCY Y( IY ) = Y( IY ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + A( L + I, J )*X( IX ) 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY 120 CONTINUE END IF END IF * RETURN * * END OF DSBMV . * END SUBROUTINE ORTHG(N,F,P,B,LDB,X,LDX,TMP) C---------------------------------------------------------------------- C C.....BLAS2 GRAM-SCHMIDT ORTHOGONALIZATION PROCEDURE FOR BLOCK LANCZOS C C THE N BY P MATRIX Z STORED IN COLUMNS F+1 TO F+P OF C ARRAY X IS REORTHOGONALIZED W.R.T. THE FIRST F COLUMNS OF C ARRAY X. THE RESULTING MATRIX Z IS THEN FACTORED INTO THE C PRODUCT OF AN N BY P ORTHONORMAL MATRIX (STORED OVER MATRIX Z) C AND A P BY P UPPER-TRIANGULAR MATRIX (STORED IN THE FIRST C P ROWS AND COLUMNS OF ARRAY B). C (BASED ON ORTHOG FROM RUTISHAUSER) C---------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,F,P,LDB,LDX REAL*8 B(LDB,1), X(LDX,1),TMP(1) INTEGER FP1,FPP LOGICAL ORIG C IF(P.EQ.0) RETURN FP1=F+1 FPP=F+P C DO 50 K=FP1,FPP ORIG=.TRUE. KM1=K-1 C C----------------------------------------------------------- 10 T=0.0D0 IF(KM1.LT.1) GO TO 25 IF(KM1.GT.1) THEN CALL DGEMV('T',N,KM1,1.0D0,X,LDX,X(1,K),1,0.0D0,TMP,1) T=T+DDOT(KM1,TMP,1,TMP,1) ELSE TMP(1)=DDOT(N,X(1,1),1,X(1,K),1) T=T+TMP(1)*TMP(1) ENDIF IF (ORIG.AND.KM1.GT.F) * CALL DCOPY(KM1-F,TMP(FP1),1,B(1,K-FP1+1),1) IF(KM1.GT.1) THEN CALL DGEMV('N',N,KM1,-1.0D0,X,LDX,TMP,1,1.0D0,X(1,K),1) ELSE CALL DAXPY(N,-TMP(1),X(1,1),1,X(1,K),1) ENDIF IF(P.EQ.1) GO TO 50 C----------------------------------------------------------- 25 S=DDOT(N,X(1,K),1,X(1,K),1) T=T+S IF(S.GT.T/100.0D0) GO TO 40 ORIG=.FALSE. GO TO 10 40 S=DSQRT(S) B(K-FP1+1,K-FP1+1)=S IF(S.NE.0.0D0) S=1.0D0/S CALL DSCAL(N,S,X(1,K),1) 50 CONTINUE RETURN END SUBROUTINE TQL2(NM,N,D,E,Z,IERR) C INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR REAL*8 D(N),E(N),Z(NM,N) REAL*8 C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2,TST1,TST2,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2, C NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND C WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD. C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS C FULL MATRIX TO TRIDIAGONAL FORM. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN C THE IDENTITY MATRIX. C C ON OUTPUT C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT C UNORDERED FOR INDICES 1,2,...,IERR-1. C C E HAS BEEN DESTROYED. C C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED C EIGENVALUES. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 TST1 = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = DABS(D(L)) + DABS(E(L)) IF (TST1 .LT. H) TST1 = H C .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... DO 110 M = L, N TST2 = TST1 + DABS(E(M)) IF (TST2 .EQ. TST1) GO TO 120 C .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP .......... 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C .......... FORM SHIFT .......... L1 = L + 1 L2 = L1 + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = PYTHAG(P,1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) D(L1) = E(L) * (P + DSIGN(R,P)) DL1 = D(L1) H = G - D(L) IF (L2 .GT. N) GO TO 145 C DO 140 I = L2, N 140 D(I) = D(I) - H C 145 F = F + H C .......... QL TRANSFORMATION .......... P = D(M) C = 1.0D0 C2 = C EL1 = E(L1) S = 0.0D0 MML = M - L C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... DO 200 II = 1, MML C3 = C2 C2 = C S2 = S I = M - II G = C * E(I) H = C * P R = PYTHAG(P,E(I)) E(I+1) = S * R S = E(I) / R C = P / R P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C .......... FORM VECTOR .......... DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C P = -S * S2 * C3 * EL1 * E(L) / DL1 E(L) = S * P D(L) = C * P TST2 = TST1 + DABS(E(L)) IF (TST2 .GT. TST1) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C .......... ORDER EIGENVALUES AND EIGENVECTORS .......... DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .GE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C .......... SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS .......... 1000 IERR = L 1001 RETURN END SUBROUTINE OPA (X,Y) C C.....MULTIPLICATION OF NROW BY NCOL SPARSE MATRIX A BY THE VECTOR X, STORE IN Y C C----------------------------------------------------------- C USER-SUPPLIED PARAMETERS FOR BLKLAN: C C SET NRMAX = ROW DIMENSION OF MATRIX A C SET NCMAX = (COLUMN DIMENSION OF MATRIX A)+1 C (ASSUMING NO. OF ROWS > NO. OF COLS) C NZMAX = MAXIMUM NO. OF NONZEROS IN MATRIX A C PARAMETER(NRMAX=500,NCMAX=200,NZMAX=2000) REAL*8 VALUE(NZMAX), ZTEMP(NRMAX,NCMAX) INTEGER POINTR(NCMAX), ROWIND(NZMAX), NROW, NCOL COMMON /MATRIXA/ NROW,NCOL,VALUE,POINTR,ROWIND,ZTEMP C----------------------------------------------------------- C COMMON /COUNT1/ MXVCOUNT C INTEGER MXVCOUNT REAL*8 X(1), Y(1) C C--------------------- C MXVCOUNT=MXVCOUNT+1 C DO 55 I=1,NROW Y(I)=0.0D0 55 CONTINUE C DO 15 I=1,NCOL C DO 10 J=POINTR(I),POINTR(I+1)-1 Y(ROWIND(J))=Y(ROWIND(J))+VALUE(J)*X(I) 10 CONTINUE C 15 CONTINUE RETURN END SUBROUTINE OPM(N,NC,X,LDX,Y,LDY) C C.....MULTIPLICATION OF A'A BY A BLOCK OF VECTORS STORED IN ARRAY X C (OUTPUT IS WRITTEN TO ARRAY Y) C C A IS AN NROW BY NCOL SPARSE MATRIX C C----------------------------------------------------------- C USER-SUPPLIED PARAMETERS FOR BLKLAN: C C SET NRMAX = ROW DIMENSION OF MATRIX A C SET NCMAX = (COLUMN DIMENSION OF MATRIX A)+1 C (ASSUMING NO. OF ROWS > NO. OF COLS) C NZMAX = MAXIMUM NO. OF NONZEROS IN MATRIX A C PARAMETER(NRMAX=500,NCMAX=200,NZMAX=2000) REAL*8 VALUE(NZMAX), ZTEMP(NRMAX,NCMAX) INTEGER POINTR(NCMAX), ROWIND(NZMAX), NROW, NCOL COMMON /MATRIXA/ NROW,NCOL,VALUE,POINTR,ROWIND,ZTEMP C----------------------------------------------------------- COMMON /COUNT1/ MXVCOUNT COMMON /COUNT2/ MTXVCOUNT C REAL*8 X(LDX,1), Y(LDY,1) C C--------------------- C DO 55 J=1,NC DO 45 I=1,NCOL Y(I,J)=0.0D0 45 CONTINUE 55 CONTINUE C C--------------------- C DO 185 J=1,NC DO 85 I=1,NROW ZTEMP(I,J)=0.0D0 85 CONTINUE 185 CONTINUE C C.....MULTIPLY BY SPARSE MATRIX A C DO 15 I=1,NCOL C DO 10 J=POINTR(I),POINTR(I+1)-1 DO 110 K=1,NC ZTEMP(ROWIND(J),K)=ZTEMP(ROWIND(J),K)+VALUE(J)*X(I,K) 110 CONTINUE 10 CONTINUE C 15 CONTINUE C C.....MULTIPLY BY TRANSPOSE (A') C DO 130 K=1,NC DO 35 I =1,NCOL DO 30 J=POINTR(I),POINTR(I+1)-1 Y(I,K)=Y(I,K)+VALUE(J)*ZTEMP(ROWIND(J),K) 30 CONTINUE 35 CONTINUE 130 CONTINUE C MXVCOUNT = MXVCOUNT+NC MTXVCOUNT=MTXVCOUNT+NC C RETURN END SUBROUTINE OPN(N, X, Y) C C.....MULTIPLICATION OF MATRIX A'A BY A VECTOR X AND STORE IN Y, WHERE C C NROW = NO. OF ROWS OF A, NCOL = NO. OF COLS OF A C C N=NCOL (ASSUME NROW << NCOL, IE. HORIZONTALLY RECTANGULAR) C C NOTE: X AND Y SHOULD BE ARRAYS OF LENGTH N C C----------------------------------------------------------- C USER-SUPPLIED PARAMETERS FOR BLKLAN: C C SET NRMAX = ROW DIMENSION OF MATRIX A C SET NCMAX = (COLUMN DIMENSION OF MATRIX A)+1 C (ASSUMING NO. OF ROWS > NO. OF COLS) C NZMAX = MAXIMUM NO. OF NONZEROS IN MATRIX A C PARAMETER(NRMAX=500,NCMAX=200,NZMAX=2000) REAL*8 VALUE(NZMAX), ZTEMP(NRMAX,NCMAX) INTEGER POINTR(NCMAX), ROWIND(NZMAX), NROW, NCOL COMMON /MATRIXA/ NROW,NCOL,VALUE,POINTR,ROWIND,ZTEMP C----------------------------------------------------------- C COMMON /COUNT1/ MXVCOUNT COMMON /COUNT2/ MTXVCOUNT C C----------------------------------------------------------- C REAL*8 X(1), Y(1) INTEGER MXVCOUNT,MTXVCOUNT C C--------------------- C DO 55 I=1,NCOL Y(I)=0.0D0 55 CONTINUE DO 85 I=1,NROW ZTEMP(I,1)=0.0D0 85 CONTINUE C C.....MULTIPLY BY SPARSE A C DO 15 I=1,NCOL C DO 10 J=POINTR(I),POINTR(I+1)-1 ZTEMP(ROWIND(J),1)=ZTEMP(ROWIND(J),1)+VALUE(J)*X(I) 10 CONTINUE C 15 CONTINUE C C.....MULTIPLY BY SPARSE A' C DO 25 I =1,NCOL C DO 20 J=POINTR(I),POINTR(I+1)-1 Y(I)=Y(I)+VALUE(J)*ZTEMP(ROWIND(J),1) 20 CONTINUE 25 CONTINUE C MXVCOUNT = MXVCOUNT+1 MTXVCOUNT=MTXVCOUNT+1 C RETURN END C SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. SCALAR ARGUMENTS .. CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC REAL*8 ALPHA, BETA * .. ARRAY ARGUMENTS .. REAL*8 A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * PURPOSE * ======= * * DGEMM PERFORMS ONE OF THE MATRIX-MATRIX OPERATIONS * * C := ALPHA*OP( A )*OP( B ) + BETA*C, * * WHERE OP( X ) IS ONE OF * * OP( X ) = X OR OP( X ) = X', * * ALPHA AND BETA ARE SCALARS, AND A, B AND C ARE MATRICES, WITH OP( A ) * AN M BY K MATRIX, OP( B ) A K BY N MATRIX AND C AN M BY N MATRIX. * * PARAMETERS * ========== * * TRANSA - CHARACTER*1. * ON ENTRY, TRANSA SPECIFIES THE FORM OF OP( A ) TO BE USED IN * THE MATRIX MULTIPLICATION AS FOLLOWS: * * TRANSA = 'N' OR 'N', OP( A ) = A. * * TRANSA = 'T' OR 'T', OP( A ) = A'. * * TRANSA = 'C' OR 'C', OP( A ) = A'. * * UNCHANGED ON EXIT. * * TRANSB - CHARACTER*1. * ON ENTRY, TRANSB SPECIFIES THE FORM OF OP( B ) TO BE USED IN * THE MATRIX MULTIPLICATION AS FOLLOWS: * * TRANSB = 'N' OR 'N', OP( B ) = B. * * TRANSB = 'T' OR 'T', OP( B ) = B'. * * TRANSB = 'C' OR 'C', OP( B ) = B'. * * UNCHANGED ON EXIT. * * M - INTEGER. * ON ENTRY, M SPECIFIES THE NUMBER OF ROWS OF THE MATRIX * OP( A ) AND OF THE MATRIX C. M MUST BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE NUMBER OF COLUMNS OF THE MATRIX * OP( B ) AND THE NUMBER OF COLUMNS OF THE MATRIX C. N MUST BE * AT LEAST ZERO. * UNCHANGED ON EXIT. * * K - INTEGER. * ON ENTRY, K SPECIFIES THE NUMBER OF COLUMNS OF THE MATRIX * OP( A ) AND THE NUMBER OF ROWS OF THE MATRIX OP( B ). K MUST * BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * ALPHA - REAL*8. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * A - REAL*8 ARRAY OF DIMENSION ( LDA, KA ), WHERE KA IS * K WHEN TRANSA = 'N' OR 'N', AND IS M OTHERWISE. * BEFORE ENTRY WITH TRANSA = 'N' OR 'N', THE LEADING M BY K * PART OF THE ARRAY A MUST CONTAIN THE MATRIX A, OTHERWISE * THE LEADING K BY M PART OF THE ARRAY A MUST CONTAIN THE * MATRIX A. * UNCHANGED ON EXIT. * * LDA - INTEGER. * ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED * IN THE CALLING (SUB) PROGRAM. WHEN TRANSA = 'N' OR 'N' THEN * LDA MUST BE AT LEAST MAX( 1, M ), OTHERWISE LDA MUST BE AT * LEAST MAX( 1, K ). * UNCHANGED ON EXIT. * * B - REAL*8 ARRAY OF DIMENSION ( LDB, KB ), WHERE KB IS * N WHEN TRANSB = 'N' OR 'N', AND IS K OTHERWISE. * BEFORE ENTRY WITH TRANSB = 'N' OR 'N', THE LEADING K BY N * PART OF THE ARRAY B MUST CONTAIN THE MATRIX B, OTHERWISE * THE LEADING N BY K PART OF THE ARRAY B MUST CONTAIN THE * MATRIX B. * UNCHANGED ON EXIT. * * LDB - INTEGER. * ON ENTRY, LDB SPECIFIES THE FIRST DIMENSION OF B AS DECLARED * IN THE CALLING (SUB) PROGRAM. WHEN TRANSB = 'N' OR 'N' THEN * LDB MUST BE AT LEAST MAX( 1, K ), OTHERWISE LDB MUST BE AT * LEAST MAX( 1, N ). * UNCHANGED ON EXIT. * * BETA - REAL*8. * ON ENTRY, BETA SPECIFIES THE SCALAR BETA. WHEN BETA IS * SUPPLIED AS ZERO THEN C NEED NOT BE SET ON INPUT. * UNCHANGED ON EXIT. * * C - REAL*8 ARRAY OF DIMENSION ( LDC, N ). * BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY C MUST * CONTAIN THE MATRIX C, EXCEPT WHEN BETA IS ZERO, IN WHICH * CASE C NEED NOT BE SET ON ENTRY. * ON EXIT, THE ARRAY C IS OVERWRITTEN BY THE M BY N MATRIX * ( ALPHA*OP( A )*OP( B ) + BETA*C ). * * LDC - INTEGER. * ON ENTRY, LDC SPECIFIES THE FIRST DIMENSION OF C AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDC MUST BE AT LEAST * MAX( 1, M ). * UNCHANGED ON EXIT. * * * LEVEL 3 BLAS ROUTINE. * * -- WRITTEN ON 8-FEBRUARY-1989. * JACK DONGARRA, ARGONNE NATIONAL LABORATORY. * IAIN DUFF, AERE HARWELL. * JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD. * SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD. * * * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. LOCAL SCALARS .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB REAL*8 TEMP * .. PARAMETERS .. REAL*8 ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. EXECUTABLE STATEMENTS .. * * SET NOTA AND NOTB AS TRUE IF A AND B RESPECTIVELY ARE NOT * TRANSPOSED AND SET NROWA, NCOLA AND NROWB AS THE NUMBER OF ROWS * AND COLUMNS OF A AND THE NUMBER OF ROWS OF B RESPECTIVELY. * NOTA = 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 IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE * * FORM C := ALPHA*A'*B + BETA*C * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN * * FORM C := ALPHA*A*B' + BETA*C * DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE * * FORM C := ALPHA*A'*B' + BETA*C * DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF * RETURN * * END OF DGEMM . * END SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. SCALAR ARGUMENTS .. REAL*8 ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. ARRAY ARGUMENTS .. REAL*8 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 - REAL*8. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * A - REAL*8 ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY A MUST * CONTAIN THE MATRIX OF COEFFICIENTS. * 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 - REAL*8 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 - REAL*8. * 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 - REAL*8 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 .. REAL*8 ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. REAL*8 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 * .. * .. 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 LOGICAL FUNCTION LSAME( CA, CB ) * * -- LAPACK AUXILIARY ROUTINE -- * ARGONNE NATIONAL LABORATORY * OCTOBER 11, 1988 * * .. SCALAR ARGUMENTS .. CHARACTER CA, CB * .. * * PURPOSE * ======= * * LSAME RETURNS .TRUE. IF CA IS THE SAME LETTER AS CB REGARDLESS * OF CASE. * * N.B. THIS VERSION OF THE ROUTINE IS ONLY CORRECT FOR ASCII CODE. * INSTALLERS MUST MODIFY THE ROUTINE FOR OTHER CHARACTER-CODES. * * FOR EBCDIC SYSTEMS THE CONSTANT IOFF MUST BE CHANGED TO -64. * FOR CDC SYSTEMS USING 6-12 BIT REPRESENTATIONS, THE SYSTEM- * SPECIFIC CODE IN COMMENTS MUST BE ACTIVATED. * * ARGUMENTS * ========= * * CA - CHARACTER*1 * CB - CHARACTER*1 * ON ENTRY, CA AND CB SPECIFY CHARACTERS TO BE COMPARED. * UNCHANGED ON EXIT. * * AUXILIARY ROUTINE FOR LEVEL 2 BLAS. * * -- WRITTEN ON 20-JULY-1986 * RICHARD HANSON, SANDIA NATIONAL LABS. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * * * .. PARAMETERS .. INTEGER IOFF PARAMETER ( IOFF = 32 ) * .. * .. INTRINSIC FUNCTIONS .. INTRINSIC ICHAR * .. * .. EXECUTABLE STATEMENTS .. * * TEST IF THE CHARACTERS ARE EQUAL * LSAME = CA.EQ.CB * * NOW TEST FOR EQUIVALENCE * IF( .NOT.LSAME ) THEN LSAME = ICHAR( CA ) - IOFF.EQ.ICHAR( CB ) END IF IF( .NOT.LSAME ) THEN LSAME = ICHAR( CA ).EQ.ICHAR( CB ) - IOFF END IF * RETURN * * THE FOLLOWING COMMENTS CONTAIN CODE FOR CDC SYSTEMS USING 6-12 BIT * REPRESENTATIONS. * * .. PARAMETERS .. * INTEGER ICIRFX * PARAMETER ( ICIRFX=62 ) * .. SCALAR ARGUMENTS .. * CHARACTER*1 CB * .. ARRAY ARGUMENTS .. * CHARACTER*1 CA(*) * .. LOCAL SCALARS .. * INTEGER IVAL * .. INTRINSIC FUNCTIONS .. * INTRINSIC ICHAR, CHAR * .. EXECUTABLE STATEMENTS .. * * SEE IF THE FIRST CHARACTER IN STRING CA EQUALS STRING CB. * * LSAME = CA(1) .EQ. CB .AND. CA(1) .NE. CHAR(ICIRFX) * * IF (LSAME) RETURN * * THE CHARACTERS ARE NOT IDENTICAL. NOW CHECK THEM FOR EQUIVALENCE. * LOOK FOR THE 'ESCAPE' CHARACTER, CIRCUMFLEX, FOLLOWED BY THE * LETTER. * * IVAL = ICHAR(CA(2)) * IF (IVAL.GE.ICHAR('A') .AND. IVAL.LE.ICHAR('Z')) THEN * LSAME = CA(1) .EQ. CHAR(ICIRFX) .AND. CA(2) .EQ. CB * END IF * * RETURN * * END OF LSAME. * END SUBROUTINE XERBLA( SRNAME, INFO ) * * -- LAPACK AUXILIARY ROUTINE -- * ARGONNE NATIONAL LABORATORY * NOVEMBER 16, 1988 * * .. 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. * * PARAMETERS * ========== * * SRNAME - CHARACTER*6. * ON ENTRY, SRNAME SPECIFIES THE NAME OF THE ROUTINE WHICH * CALLED XERBLA. * * INFO - INTEGER. * ON ENTRY, INFO SPECIFIES 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 C C @(#)RANDOM.F 3.2 (BNP) 12/9/88; FROM RANDOM.F 2.2 10/13/87 C REAL*8 FUNCTION RANDOM(IY) INTEGER IY C C RANDOM IS A UNIFORM RANDOM NUMBER GENERATOR BASED ON THEORY AND C SUGGESTIONS GIVEN IN D.E.KNUTH (1969), VOL 2. THE INTEGER IY C SHOULD BE INITIALIZED TO AN ARBITRARY INTEGER PRIOR TO THE FIRST CALL C TO RANDOM.THE CALLING PROGRAM SHOULD NOT ALTER THE VALUE OF IY C BETWEEN SUBSEQUENT CALLS TO RANDOM.VALUES OF RANDOM WILL BE RETURNED C IN THE INTERVAL (0,1). C INTEGER M,IA,IC,MIC REAL*8 HALFM,S C INTEGER M2,ITWO DATA M2,ITWO/0,2/ IF (M2.EQ.0) THEN C C IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH C M = 1 10 M2 = M M = ITWO*M2 IF (M.GT.M2) GO TO 10 HALFM = M2 C C COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD C IA = 8*IDINT(HALFM*DATAN(1.D0)/8.D0)+5 IC = 2*IDINT(HALFM*(0.5D0-DSQRT(3.D0)/6.D0))+1 MIC = (M2-IC)+M2 C C S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT C S = 0.5D0/HALFM C C COMPUTE NEXT RANDOM NUMBER C ENDIF IY = IY*IA C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHICH DO NOT ALLOW C INTEGER OVERFLOW ON ADDITION C IF (IY.GT.MIC) IY = (IY-M2)-M2 C IY = IY+IC C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE THE C WORD LENGTH FOR ADDITION IS GREATER THAN FOR MULTIPLICATION C IF (IY/2.GT.M2) IY = (IY-M2)-M2 C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER C OVERFLOW AFFECTS THE SIGN BIT C IF (IY.LT.0) IY = (IY+M2)+M2 RANDOM = DBLE(IY)*S RETURN END REAL*8 FUNCTION PYTHAG(A,B) REAL*8 A,B C C FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW C REAL*8 P,R,S,T,U C C DECLARE REVISION CONTROL SYSTEM (RCS) HEADER C CHARACTER*80 HDR DATA HDR /'$HEADER: PYTHAG.F,V 1.1 86/03/03 08:52:24 BERRY EXP $'/ C P = DMAX1(DABS(A),DABS(B)) IF (P .EQ. 0.0D0) GO TO 20 R = (DMIN1(DABS(A),DABS(B))/P)**2 10 CONTINUE T = 4.0D0 + R IF (T .EQ. 4.0D0) GO TO 20 S = R/T U = 1.0D0 + 2.0D0*S P = U*P R = (S/U)**2 * R GO TO 10 20 PYTHAG = P RETURN END .