C************************************************************************* C * C SPARSE SVD VIA HYBRID BLOCK LANCZOS METHOD ON 2-CYCLIC 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 BLS1 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) REAL*8 TMPU(LDU),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 NROW, NCOL, 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 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(NCMAX=500,NZMAX=2000) REAL*8 VALUE(NZMAX) INTEGER POINTR(NCMAX), ROWIND(NZMAX) COMMON /MATRIXA/ VALUE,POINTR,ROWIND C----------------------------------------------------------- C MATRIX=7 IN=5 RESULTS=9 IVEC=8 MXVCOUNT=0 MTXVCOUNT=0 C C**************************************************************** C C OPEN FILES FOR INPUT/OUTPUT C OPEN(IN, FILE='BLI5') OPEN(MATRIX, FILE='LAI7') OPEN(RESULTS,FILE='BLO9') C READ(IN,*) NAME, NC, NB, NUMS, TOL, VECTORS IF(VECTORS) OPEN(IVEC,FILE='BLO8',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.GE. NCMAX .OR. NNZERO .GT. NZMAX) 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(M,N,NKOLD,LDV,V,LDU,U,S,NCOLD,NBOLD,ITER, * TOL,RES,150,NK,NC,NB,IMEM) C TIME=ETIME(TMP2(1))-T1 C C.......WRITE UNFORMATTED FILE CONTAINING LEFT AND RIGHT S-VECTOR PAIRS C IF(VECTORS) THEN DO 510 J = 1,NK TMP1=DDOT(N,V(1,J),1,V(1,J),1) TMP3=DDOT(M,U(1,J),1,U(1,J),1) TMP4=DSQRT(TMP1+TMP3) CALL OP(M,N,V(1,J),TMPU) CALL OPT(N,M,U(1,J),TMPV) CALL DAXPY(M,-S(J),U(1,J),1,TMPU,1) CALL DAXPY(N,-S(J),V(1,J),1,TMPV,1) XNORM=DDOT(N,TMPV,1,TMPV,1)+ * DDOT(M,TMPU,1,TMPU,1) XNORM=DSQRT(XNORM) RES(J)=XNORM/TMP4 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 (CYCLIC)' * /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 C C******************************************************************* C SUBROUTINE BLKLAN(M,N,IK,LDV,V,LDU,U, * SING,IC,IB,ITER,TOL,RES,MAXIT, * IKO,ICO,IBO,IMEM) C IMPLICIT REAL*8 (A-H,O-Z) INTEGER M,N,LDV,LDU,IB,IC,IK,ITER,MAXIT REAL*8 V(LDV,IK),U(LDU,IK),SING(IK),RES(IK) REAL*8 RANDOM C C *************************************** C * * C * [BLOCK LANCZOS SVD ALGORITHM] * C * * C * FIVE PHASES: * C * * C * PHASE 1: BLOCK LANCZOS TO YIELD * C * BLOCK UPPER-TRI. MATRIX S WHICH * C * SHARES THE SAME S-VALUES AS A. * C * TOTAL RE-ORTHOGONALIZATION USED * C * * C * PHASE 2: LANCZOS METHOD FOR BI- * C * DIAGONALIZING THE S MATRIX FROM * C * PHASE 1 TO YIELD THE BI-DIAGONAL * C * MATRIX B WHICH PRESERVES THE * C * SAME S-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 QR * C * ITERATION TO DIAGONALIZE B AND * C * HENCE PRODUCE THE APPROX. S-VALUES * C * (ARRAY ALPHA) OF MATRIX A. * C * * C * PHASE 4: CONVERGENCE TEST. * C * * C * PHASE 5: ITERATION RESTART. * C * PROJECT AGAINST CONVERGED * C * S-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 TRIPLETS * C * TOL = RESIDUAL TOLERANCE * C * MAXIT = MAXIMUM # ITERATIONS * C * * C * OUTPUT: ITER = 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. TRIPLETS FOUND * C * SING(IKO) = S-VALUES * C * U(LDU,IKO) = LEFT S-VECTORS * C * V(LDV,IKO) = RIGHT S-VECTORS * C * * C * REQUIRED: * C * * C * CALLS OP(M,N,X,Y) FOR Y=A*X * C * CALLS OPT(N,M,X,Y) FOR Y=A'*X * C * * C * EXTERNALS ROUTINES: * C * * C * RANDOM - RANDOM NO. GENERATOR * C * BLAS1 KERNELS USED: * C * DNRM2 - EUCLIDEAN VEC. NORM * C * DSCAL - VECTOR SCALING * C * DSWAP - VECTOR SWAP * C * DCOPY - VECTOR COPY * C * DDOT - DOT PRODUCT * C * DAXPY - VECTOR OPERATION * C * DROT - PLANE ROTATIONS * C * DROTG - GIVENS ROTATIONS * C * * C * BLAS2 AND 3 KERNELS USED: * C * DGEMM - C=C-A*B (LAPACK BLAS3) * C * DGEMV - Y=Y-A*X (LAPACK BLAS2) * C * ORTHG - MOD. GRAM-SCHMIDT * C * DTBMV - BAND MATRIX-VECTOR MULT.* C * (LAPACK BLAS2) * C * * C *************************************** C C----------------------------------------------------------- C USER-SUPPLIED PARAMETERS FOR BLKLAN: C 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(NCMAX=500,NZMAX=2000) REAL*8 VALUE(NZMAX) INTEGER POINTR(NCMAX), ROWIND(NZMAX) COMMON /MATRIXA/ VALUE,POINTR,ROWIND C 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 M2 = ROW DIMENSION OF MATRIX A (M) C N2 = COL DIMENSION OF MATRIX A (N) C LDT2 = MAX(M2,N2) C--------------------------------------------------------- C PARAMETER( * IB2=10, * IK2=10, * IC2=50, * M2=374, * N2=82, * LDT2=374) C REAL*8 W(M2,IB2),Y(N2,IB2),S(IB2,IC2),R(IB2,IC2-IB2), * TEMP(LDT2,IB2),UU(M2,IC2),VV(N2,IC2),U0(M2,IK2), * V0(N2,IK2),TRES(IB2),UVTMP(LDT2,IC2+IK2), * BIGS(IB2+1,IC2),P(N2),Q(M2),T(M2),PP(IC2,IC2), * QQ(IC2,IC2),ALPHA(IC2),Z(N2),BETA(IC2) C INTEGER KINDX(IB2) C C--------------------------------------------------------- C C NS IS THE NUMBER OF BLOCK COLUMN PARTITIONS TO FORM C C NB IS THE NUMBER OF COLUMNS IN U(J), V(J) C C I.E., V = [ V(1), V(2), ... , V(NS) ], AND C U = [ U(1), U(2), ... , U(NS) ]. C C......SET DIMENSIONAL PARAMETERS (LOCAL) C NB=IB NC=IC K=IK C C----------------------------------------------------------------------- C C......REAL STORAGE REQUIRED C IMEM=LDV*IK+LDU*IK+2*IK+ * M2*IB2+N2*IB2+IB2*IC2+IB2*(IC2-IB2)+ * M2*IK2+N2*IK2+LDT2*IB2+ * M2*IC2+N2*IC2+ * LDT2*(IC2+IK2)+ * (IB2+1)*IC2+N2+2*M2+IC2-1+ * 2*IC2*IC2+IC2+IB2+N2+IC2-1+2*IC2 IMEM=8*IMEM+4*IB2 C C.....IMEM STORES TOTAL AUXILIARY STORAGE REQUIRED (IN BYTES) 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 C ******************** B E G I N P H A S E 1 ********************** C C......DETERMINE NUMBER OF BLOCKS FOR THIS 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......INITIALIZE S AND R ARRAYS 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 W(1) = A*V(1) C DO 5 J=1,NB CALL OP(M,N,VV(1,J),W(1,J)) 5 CONTINUE C IF(ICONV.GT.0) THEN C C...............ORTHOGONALIZE W(1) W.R.T. U(0) C DO 675 JJ=1,ICONV DO 676 II=1,M UVTMP(II,JJ)=U0(II,JJ) 676 CONTINUE 675 CONTINUE C ENDIF C DO 775 JJ=1,NB DO 776 II=1,M UVTMP(II,ICONV+JJ)=W(II,JJ) W(II,JJ)=0.0D0 776 CONTINUE 775 CONTINUE C CALL ORTHG(M,ICONV,NB,W,M2,UVTMP,LDT2,TEMP) C C LOAD S(J) SUBMATRIX C DO 160 JJ = 1,NB DO 170 I= 1,JJ S(I,JJ)=W(I,JJ) 170 CONTINUE 160 CONTINUE C C LOAD U(J) SUBMATRIX C DO 512 JJ=1,NB DO 513 II=1,M UU(II,JJ)=UVTMP(II,JJ+ICONV) 513 CONTINUE 512 CONTINUE 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^T * U(J) - V(J) * S(J)^T C DO 64 JJ=1,NB CALL OPT(N,M,UU(1,JPTR+JJ-1),Y(1,JJ)) 64 CONTINUE C CALL DGEMM('N','T',N,NB,NB,-1.0D0,VV(1,JPTR),N2, * S(1,JPTR),LDS,1.0D0,Y,N2) 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 IF(KK.LE.NB*(J-1)) THEN DO 732 LL=1,N UVTMP(LL,KK)=VV(LL,KK) 732 CONTINUE ELSE DO 733 LL=1,N UVTMP(LL,KK)=V0(LL,KK-NB*(J-1)) 733 CONTINUE ENDIF 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 C 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 CALL ORTHG(N,NJ,NB,Y,N2,UVTMP,LDT2,TEMP) C C GET QR FACTORIZATION: Y(J) = V(J) * R(J), [ R(J) := UPPER TRI.] C C INCREMENT POINTER FOR V(J) C JINC=JPTR+NB 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 V(J) SUBMATRIX C DO 610 JJ=1,NB DO 611 II=1,N VV(II,JINC+JJ-1)=UVTMP(II,NJ+JJ) 611 CONTINUE 610 CONTINUE C C COMPUTE W(J) = A * V(J) - U(J) * R(J)^T C DO 105 JJ=1,NB CALL OP(M,N,VV(1,JINC+JJ-1),W(1,JJ)) 105 CONTINUE C CALL DGEMM('N','T',M,NB,NB,-1.0D0,UU(1,JPTR),M2, * R(1,JPTR),LDR,1.0D0,W,M2) C C-------------------------------------------------------------- C..........ORTHOGONALIZE W(J) W.R.T. [U(1),....U(J) | U(0) ] C-------------------------------------------------------------- C IF(ICONV.GT.0) THEN C C...............ORTHOGONALIZE W(J) W.R.T. [U(0), U(1), ... , U(J)] C NJ=NB*(J-1)+ICONV DO 831 KK=1,NJ IF(KK.LE.NB*(J-1)) THEN DO 832 LL=1,M UVTMP(LL,KK)=UU(LL,KK) 832 CONTINUE ELSE DO 835 LL=1,M UVTMP(LL,KK)=U0(LL,KK-NB*(J-1)) 835 CONTINUE ENDIF 831 CONTINUE C ELSE C C...............ORTHOGONALIZE W(J+1) W.R.T. [U(1), ... ,U(J)] C NJ=NB*(J-1) DO 833 KK=1,NJ DO 834 LL=1,M UVTMP(LL,KK)=UU(LL,KK) 834 CONTINUE 833 CONTINUE ENDIF DO 931 KK=1,NB DO 932 LL=1,M UVTMP(LL,NJ+KK)=W(LL,KK) W(LL,KK)=0.0D0 932 CONTINUE 931 CONTINUE C CALL ORTHG(M,NJ,NB,W,M2,UVTMP,LDT2,TEMP) C C GET QR FACTORIZATION: W(J) = U(J) * S(J), [ S(J) := UPPER TRI.] C C LOAD S(J) SUBMATRIX C DO 140 JJ = 1,NB DO 150 I= 1,JJ S(I,JJ+JINC-1)=W(I,JJ) 150 CONTINUE 140 CONTINUE C C LOAD U(J) SUBMATRIX C DO 710 JJ=1,NB DO 711 II=1,M UU(II,JINC+JJ-1)=UVTMP(II,NJ+JJ) 711 CONTINUE 710 CONTINUE C 30 CONTINUE C C----------------------------------------------------------------- C FORM THE BLOCK UPPER-TRIANGULAR S MATRIX IN BIGS ARRAY C----------------------------------------------------------------- C C....ARRAY PP WILL STORE LEFT S-VECTORS OF [BIGS] <-BANDED STRUCTURE C....ARRAY QQ WILL STORE RIGHT S-VECTORS OF [BIGS] 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 BI-DIAGONALIZATION: PP GENERATES RIGHT S-VECTORS OF BAND MATRIX C QQ GENERATES LEFT S-VECTORS OF BAND MATRIX C C CHOOSE A STARTING VECTOR P(1) [HAVING UNITY 2-NORM] C IHBW=NB C C DO 200 I=1,NN P(I)=RANDOM(IRAND) 200 CONTINUE C PNM=DNRM2(NN,P,1) C DO 201 I=1,NN P(I)=P(I)/PNM 201 CONTINUE C DO 205 I=1,NN T(I)=P(I) 205 CONTINUE C DO 206 II = 1,NN PP(II,1)=P(II) 206 CONTINUE C C COMPUTE T(1) := S * P(1) C CALL DTBMV('U','N','N',NN,IHBW,BIGS,LDB,T,1) C C COMPUTE ALPHA(1) C ALPHA(1)=DNRM2(NN,T,1) C C ITERATE FOR J = 1, 2, ... , NN C DO 210 J = 1,NN C IF(ALPHA(J).NE.0.0D0) THEN C C COMPUTE Q(J) := T(J)/ALPHA(J) C DO 220 JJ = 1,NN Q(JJ)=T(JJ)/ALPHA(J) Z(JJ)=Q(JJ) 220 CONTINUE DO 225 II = 1,NN QQ(II,J)=Q(II) 225 CONTINUE C IF(J.EQ.NN) GO TO 500 C C COMPUTE Z(J) := S^T * Q(J) - ALPHA(J) * P(J) C CALL DTBMV('U','T','N',NN,IHBW,BIGS,LDB,Z,1) 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 C IF(J.GT.1) THEN C DO 1626 II=1,NN UVTMP(II,J+1)=Z(II) 1626 CONTINUE CALL ORTHG(NN,J,1,W,M2,UVTMP,LDT2,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 779 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 P(JJ)=Z(JJ)/BETA(J) T(JJ)=P(JJ) 230 CONTINUE DO 235 II = 1,NN PP(II,J+1)=P(II) 235 CONTINUE C C COMPUTE T(J+1) := S * P(J+1) - BETA(J) * Q(J) C CALL DTBMV('U','N','N',NN,IHBW,BIGS,LDB,T,1) CALL DAXPY(NN,-BETA(J),Q,1,T,1) C C----------------------------------------------------------------------- C.......ORTHOGONALIZE T W.R.T. PREVIOUS QQ'S C----------------------------------------------------------------------- C DO 1720 JJ=1,J DO 1725 II=1,NN UVTMP(II,JJ)=QQ(II,JJ) 1725 CONTINUE 1720 CONTINUE C IF(J.GT.1) THEN C DO 1721 II=1,NN UVTMP(II,J+1)=T(II) 1721 CONTINUE CALL ORTHG(NN,J,1,W,M2,UVTMP,LDT2,TEMP) DO 1722 II=1,NN T(II)=UVTMP(II,J+1) 1722 CONTINUE ELSE DUM=-DDOT(NN,UVTMP(1,1),1,T,1) CALL DAXPY(NN,DUM,UVTMP(1,1),1,T,1) ENDIF C C COMPUTE ALPHA(J+1) C ALPHA(J+1)=DNRM2(NN,T,1) C ELSE C GO TO 500 C ENDIF C ELSE C GO TO 500 C ENDIF C 210 CONTINUE C C ******************** B E G I N P H A S E 2A (NB=1 CASE) ********** C C (POINT ALGORITHM) C C LANCZOS BI-DIAGONALIZATION: VV GENERATES RIGHT S-VECTORS (VL=N) C UU GENERATES LEFT S-VECTORS (VL=M) C COMPUTE T(1) := A * P(1) C 888 CALL OP(M,N,P,T) C C----------------------------------------------------------------- C.....ORTHOGONALIZE T W.R.T. CONVERGED LEFT S-VECTORS (U0 ARRAY) C----------------------------------------------------------------- C IF(ICONV.GT.0) THEN C DO 600 JJ=1,ICONV DO 605 II=1,M UVTMP(II,JJ)=U0(II,JJ) 605 CONTINUE 600 CONTINUE C IF(ICONV.GT.1) THEN DO 1821 II=1,M UVTMP(II,ICONV+1)=T(II) 1821 CONTINUE CALL ORTHG(M,ICONV,1,W,M2,UVTMP,LDT2,TEMP) DO 1822 II=1,M T(II)=UVTMP(II,ICONV+1) 1822 CONTINUE ELSE DUM=-DDOT(NN,UVTMP(1,1),1,T,1) CALL DAXPY(NN,DUM,UVTMP(1,1),1,T,1) ENDIF ENDIF C C----------------------------------------------------------------- C C COMPUTE ALPHA(1) C ALPHA(1)=DNRM2(M,T,1) C C ITERATE FOR J =1, ... , NN C DO 310 J = 1,NN C IF(ALPHA(J).NE.0.0D0) THEN C C COMPUTE Q(J) := T(J)/ALPHA(J) C DO 320 JJ = 1,M Q(JJ)=T(JJ)/ALPHA(J) 320 CONTINUE DO 321 JJ = 1,M UU(JJ,J)=Q(JJ) 321 CONTINUE C NVECU=J IF(J.EQ.NN) GOTO 500 C C COMPUTE Z(J) := A^T * Q(J) - ALPHA(J) * P(J) C CALL OPT(N,M,Q,Z) 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 326 II=1,M U0(II,ICONV+1)=Q(II) 326 CONTINUE DO 327 II=1,N V0(II,ICONV+1)=P(II) 327 CONTINUE SING(ICONV+1)=ALPHA(1) RES(ICONV+1)=ZNM ICONV=ICONV+1 K=K-1 IF(K.EQ.0) THEN ITER=ITER-1 GO TO 999 ENDIF 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 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,M2,UVTMP,LDT2,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 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 P(JJ)=Z(JJ)/BETA(J) 330 CONTINUE DO 331 JJ = 1,N VV(JJ,J+1)=P(JJ) 331 CONTINUE NVECV=J+1 C C COMPUTE T(J+1) := A * P(J+1) - BETA(J) * Q(J) C CALL OP(M,N,P,T) CALL DAXPY(M,-BETA(J),Q,1,T,1) C C----------------------------------------------------------------------- C.......ORTHOGONALIZE T W.R.T. CONVERGED LEFT S-VECTORS AND PREVIOUS UU'S C----------------------------------------------------------------------- C DO 720 JJ=1,ICONV+J DO 725 II=1,M IF(JJ.LE.ICONV) THEN UVTMP(II,JJ)=U0(II,JJ) ELSE UVTMP(II,JJ)=UU(II,JJ-ICONV) ENDIF 725 CONTINUE 720 CONTINUE C IF(ICONV+J.GT.1) THEN DO 2021 II=1,M UVTMP(II,ICONV+J+1)=T(II) 2021 CONTINUE CALL ORTHG(M,ICONV+J,1,W,M2,UVTMP,LDT2,TEMP) DO 2022 II=1,M T(II)=UVTMP(II,ICONV+J+1) 2022 CONTINUE ELSE DUM=-DDOT(M,UVTMP(1,1),1,T,1) CALL DAXPY(M,DUM,UVTMP(1,1),1,T,1) ENDIF C C COMPUTE ALPHA(J+1) C ALPHA(J+1)=DNRM2(M,T,1) C ELSE C GO TO 500 C ENDIF C ELSE C GO TO 500 C ENDIF C 310 CONTINUE C 500 CONTINUE C C ******************** B E G I N P H A S E 3 ********************** C IF(NB.GT.1) THEN CALL QRITER2(NN,NN,ALPHA,BETA,QQ,IC2,PP,IC2,INFO) ELSE C DO 735 JJ=1,NN DO 736 II=1,NN PP(II,JJ)=0.0D0 QQ(II,JJ)=0.0D0 736 CONTINUE 735 CONTINUE DO 737 II=1,NN PP(II,II)=1.0D0 QQ(II,II)=1.0D0 737 CONTINUE C CALL QRITER2(NN,NN,ALPHA,BETA,QQ,IC2,PP,IC2,INFO) C ENDIF 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 SING(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.EQ.0) GO TO 997 C C.........DETERMINE NUMBER OF BLOCKS FOR THIS 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 C CALL MXMA(UU,1,M2,QQ,1,IC2,U,1,LDU,M,NN,NB+K0) CALL DGEMM('N','N',M,NB+K0,NN,1.0D0,UU,M2,QQ,IC2, * 0.0D0,U,LDU) C CALL MXMA(VV,1,N2,PP,1,IC2,V,1,LDV,N,NN,NB+K0) CALL DGEMM('N','N',N,NB+K0,NN,1.0D0,VV,N2,PP,IC2, * 0.0D0,V,LDV) ELSE C CALL MXVA(UU,1,M2,QQ(1,1),1,U(1,1),1,M,NN) CALL DGEMV('N',M,NN,1.0D0,UU,M2,QQ(1,1),1,0.0D0, * U(1,1),1) C CALL MXVA(VV,1,N2,PP(1,1),1,V(1,1),1,N,NN) CALL DGEMV('N',N,NN,1.0D0,VV,N2,PP(1,1),1,0.0D0, * V(1,1),1) ENDIF C IF(K0.GT.0) THEN DO 420 JJ=ICONV+1,ICONV+K0 DO 425 II=1,M U0(II,JJ)=U(II,KINDX(JJ-ICONV)) 425 CONTINUE 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 RIGHT S-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 925 II=1,M U(II,JJ)=U0(II,JJ) 925 CONTINUE 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 OP (M, N, X, Y) C C.....MULTIPLICATION OF M BY N MATRIX A BY A VECTOR X, STORE IN Y C C----------------------------------------------------------- C USER-SUPPLIED PARAMETERS FOR BLKLAN: C 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(NCMAX=500,NZMAX=2000) REAL*8 VALUE(NZMAX) INTEGER POINTR(NCMAX), ROWIND(NZMAX) COMMON /MATRIXA/ VALUE,POINTR,ROWIND C C----------------------------------------------------------- C COMMON /COUNT1/ MXVCOUNT C INTEGER N,MXVCOUNT REAL*8 X(1), Y(1) C C--------------------- C MXVCOUNT=MXVCOUNT+1 C DO 55 I=1,M Y(I)=0.0D0 55 CONTINUE C DO 15 I=1,N 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 OPT (N, M, X, Y) C C.....MULTIPLICATION OF N BY M MATRIX A' BY A VECTOR X, STORE IN Y C C----------------------------------------------------------- C USER-SUPPLIED PARAMETERS FOR BLKLAN: C 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(NCMAX=500,NZMAX=2000) REAL*8 VALUE(NZMAX) INTEGER POINTR(NCMAX), ROWIND(NZMAX) COMMON /MATRIXA/ VALUE,POINTR,ROWIND C C----------------------------------------------------------- COMMON /COUNT2/ MTXVCOUNT C INTEGER N,MTXVCOUNT REAL*8 X(1), Y(1) C C--------------------- C MTXVCOUNT=MTXVCOUNT+1 C DO 55 I=1,N Y(I)=0.0D0 55 CONTINUE C DO 25 I =1,N C DO 20 J=POINTR(I),POINTR(I+1)-1 Y(I)=Y(I)+VALUE(J)*X(ROWIND(J)) 20 CONTINUE 25 CONTINUE C RETURN END C 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 QRITER2(N,P,S,E,U,LDU,V,LDV,INFO) IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,P,LDU,LDV,INFO REAL*8 S(N),E(N),U(LDU,1),V(LDV,1) C C ************************************** C * * C * REDUCTION FROM BI-DIAGONAL TO * C * DIAGONAL FORM IN PHASE 3 OF * C * ROUTINE BLKLAN2 (QR ITERATION). * C * * C * THIS IS COPIED FROM THE LINPACK * C * ROUTINE DSVDC. VECTORS * C * ARE ACCUMULATED FOR THIS VERSION.* C * * C * ON INPUT, * C * ARRAY S CONTAINS THE DIAGONAL * C * ELEMENTS OF MATRIX B AND E * C * CONTAINS THE SUP-DIAGONAL. * C * * C * ON OUTPUT, * C * ARRAY S CONTAINS THE SINGULAR * C * VALUES OF THE ORIGINAL MATRIX A. * C * THE VARIABLE ITER RETURNS THE * C * NUMBER OF QR ITERATIONS REQUIRED.* C * * C * ARRAY U:= LEFT S-VECTORS, * C * ARRAY V:= RIGHT S-VECTORS. * C * * C ************************************** C C******************************************************** C INTERNALS C LOGICAL WANTU,WANTV C C MAIN ITERATION LOOP FOR THE SINGULAR VALUES. C WANTU=.TRUE. WANTV=.TRUE. C MAXIT=30 M=MIN0(P,N+1) MM = M ITER = 0 360 CONTINUE C C QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND. C C ...EXIT IF (M .EQ. 0) GO TO 620 C C IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET C FLAG AND RETURN. C IF (ITER .LT. MAXIT) GO TO 370 INFO = M C ......EXIT GO TO 620 370 CONTINUE C C THIS SECTION OF THE PROGRAM INSPECTS FOR C NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS. ON C COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS. C C KASE = 1 IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M C KASE = 2 IF S(L) IS NEGLIGIBLE AND L.LT.M C KASE = 3 IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND C S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP). C KASE = 4 IF E(M-1) IS NEGLIGIBLE (CONVERGENCE). C DO 390 LL = 1, M L = M - LL C ...EXIT IF (L .EQ. 0) GO TO 400 TEST = DABS(S(L)) + DABS(S(L+1)) ZTEST = TEST + DABS(E(L)) IF (ZTEST .NE. TEST) GO TO 380 E(L) = 0.0D0 C ......EXIT GO TO 400 380 CONTINUE 390 CONTINUE 400 CONTINUE IF (L .NE. M - 1) GO TO 410 KASE = 4 GO TO 480 410 CONTINUE LP1 = L + 1 MP1 = M + 1 DO 430 LLS = LP1, MP1 LS = M - LLS + LP1 C ...EXIT IF (LS .EQ. L) GO TO 440 TEST = 0.0D0 IF (LS .NE. M) TEST = TEST + DABS(E(LS)) IF (LS .NE. L + 1) TEST = TEST + DABS(E(LS-1)) ZTEST = TEST + DABS(S(LS)) IF (ZTEST .NE. TEST) GO TO 420 S(LS) = 0.0D0 C ......EXIT GO TO 440 420 CONTINUE 430 CONTINUE 440 CONTINUE IF (LS .NE. L) GO TO 450 KASE = 3 GO TO 470 450 CONTINUE IF (LS .NE. M) GO TO 460 KASE = 1 GO TO 470 460 CONTINUE KASE = 2 L = LS 470 CONTINUE 480 CONTINUE L = L + 1 C C PERFORM THE TASK INDICATED BY KASE. C GO TO (490,520,540,570), KASE C C DEFLATE NEGLIGIBLE S(M). C 490 CONTINUE MM1 = M - 1 F = E(M-1) E(M-1) = 0.0D0 DO 510 KK = L, MM1 K = MM1 - KK + L T1 = S(K) CALL DROTG(T1,F,CS,SN) S(K) = T1 IF (K .EQ. L) GO TO 500 F = -SN*E(K-1) E(K-1) = CS*E(K-1) 500 CONTINUE IF (WANTV) CALL DROT(P,V(1,K),1,V(1,M),1,CS,SN) 510 CONTINUE GO TO 610 C C SPLIT AT NEGLIGIBLE S(L). C 520 CONTINUE F = E(L-1) E(L-1) = 0.0D0 DO 530 K = L, M T1 = S(K) CALL DROTG(T1,F,CS,SN) S(K) = T1 F = -SN*E(K) E(K) = CS*E(K) IF (WANTU) CALL DROT(N,U(1,K),1,U(1,L-1),1,CS,SN) 530 CONTINUE GO TO 610 C C PERFORM ONE QR STEP. C 540 CONTINUE C C CALCULATE THE SHIFT. C SCALE = DMAX1(DABS(S(M)),DABS(S(M-1)),DABS(E(M-1)), * DABS(S(L)),DABS(E(L))) SM = S(M)/SCALE SMM1 = S(M-1)/SCALE EMM1 = E(M-1)/SCALE SL = S(L)/SCALE EL = E(L)/SCALE B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0D0 C = (SM*EMM1)**2 SHIFT = 0.0D0 IF (B .EQ. 0.0D0 .AND. C .EQ. 0.0D0) GO TO 550 SHIFT = DSQRT(B**2+C) IF (B .LT. 0.0D0) SHIFT = -SHIFT SHIFT = C/(B + SHIFT) 550 CONTINUE F = (SL + SM)*(SL - SM) + SHIFT G = SL*EL C C CHASE ZEROS. C MM1 = M - 1 DO 560 K = L, MM1 CALL DROTG(F,G,CS,SN) IF (K .NE. L) E(K-1) = F F = CS*S(K) + SN*E(K) E(K) = CS*E(K) - SN*S(K) G = SN*S(K+1) S(K+1) = CS*S(K+1) IF (WANTV) CALL DROT(P,V(1,K),1,V(1,K+1),1,CS,SN) CALL DROTG(F,G,CS,SN) S(K) = F F = CS*E(K) + SN*S(K+1) S(K+1) = -SN*E(K) + CS*S(K+1) G = SN*E(K+1) E(K+1) = CS*E(K+1) IF (WANTU .AND. K .LT. N) * CALL DROT(N,U(1,K),1,U(1,K+1),1,CS,SN) 560 CONTINUE E(M-1) = F ITER = ITER + 1 GO TO 610 C C CONVERGENCE. C 570 CONTINUE C C MAKE THE SINGULAR VALUE POSITIVE. C IF (S(L) .GE. 0.0D0) GO TO 580 S(L) = -S(L) IF (WANTV) CALL DSCAL(P,-1.0D0,V(1,L),1) 580 CONTINUE C C ORDER THE SINGULAR VALUE. C 590 IF (L .EQ. MM) GO TO 600 C ...EXIT IF (S(L) .GE. S(L+1)) GO TO 600 T = S(L) S(L) = S(L+1) S(L+1) = T IF (WANTV .AND. L .LT. P) * CALL DSWAP(P,V(1,L),1,V(1,L+1),1) IF (WANTU .AND. L .LT. N) * CALL DSWAP(N,U(1,L),1,U(1,L+1),1) L = L + 1 GO TO 590 600 CONTINUE ITER = 0 M = M - 1 610 CONTINUE GO TO 360 620 CONTINUE RETURN END SUBROUTINE FORMBIGS(N,NDP,R,S,BIGS,LDB,LDR,LDS) INTEGER N,LDB,LDR,LDS,NDP 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 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 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 DTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) * .. SCALAR ARGUMENTS .. INTEGER INCX, K, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. ARRAY ARGUMENTS .. REAL*8 A( LDA, * ), X( * ) * .. * * PURPOSE * ======= * * DTBMV PERFORMS ONE OF THE MATRIX-VECTOR OPERATIONS * * X := A*X, OR X := A'*X, * * WHERE X IS AN N ELEMENT VECTOR AND A IS AN N BY N UNIT, OR NON-UNIT, * UPPER OR LOWER TRIANGULAR BAND MATRIX, WITH ( K + 1 ) DIAGONALS. * * PARAMETERS * ========== * * UPLO - CHARACTER*1. * ON ENTRY, UPLO SPECIFIES WHETHER THE MATRIX IS AN UPPER OR * LOWER TRIANGULAR MATRIX AS FOLLOWS: * * UPLO = 'U' OR 'U' A IS AN UPPER TRIANGULAR MATRIX. * * UPLO = 'L' OR 'L' A IS A LOWER TRIANGULAR MATRIX. * * UNCHANGED ON EXIT. * * TRANS - CHARACTER*1. * ON ENTRY, TRANS SPECIFIES THE OPERATION TO BE PERFORMED AS * FOLLOWS: * * TRANS = 'N' OR 'N' X := A*X. * * TRANS = 'T' OR 'T' X := A'*X. * * TRANS = 'C' OR 'C' X := A'*X. * * UNCHANGED ON EXIT. * * DIAG - CHARACTER*1. * ON ENTRY, DIAG SPECIFIES WHETHER OR NOT A IS UNIT * TRIANGULAR AS FOLLOWS: * * DIAG = 'U' OR 'U' A IS ASSUMED TO BE UNIT TRIANGULAR. * * DIAG = 'N' OR 'N' A IS NOT ASSUMED TO BE UNIT * TRIANGULAR. * * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE ORDER OF THE MATRIX A. * N MUST BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * K - INTEGER. * ON ENTRY WITH UPLO = 'U' OR 'U', K SPECIFIES THE NUMBER OF * SUPER-DIAGONALS OF THE MATRIX A. * ON ENTRY WITH UPLO = 'L' OR 'L', K SPECIFIES THE NUMBER OF * SUB-DIAGONALS OF THE MATRIX A. * K MUST SATISFY 0 .LE. K. * 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 MATRIX OF COEFFICIENTS, 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 AN UPPER * TRIANGULAR 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 MATRIX OF COEFFICIENTS, 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 A LOWER * TRIANGULAR 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 * * NOTE THAT WHEN DIAG = 'U' OR 'U' THE ELEMENTS OF THE ARRAY A * CORRESPONDING TO THE DIAGONAL ELEMENTS OF THE MATRIX ARE NOT * REFERENCED, 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 * ( 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 N * ELEMENT VECTOR X. ON EXIT, X IS OVERWRITTEN WITH THE * TRANFORMED VECTOR X. * * INCX - INTEGER. * ON ENTRY, INCX SPECIFIES THE INCREMENT FOR THE ELEMENTS OF * X. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. REAL*8 ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. REAL*8 TEMP INTEGER I, INFO, IX, J, JX, KPLUS1, KX, L LOGICAL NOUNIT * .. 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( .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( K.LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.( K + 1 ) )THEN INFO = 7 ELSE IF( INCX.EQ.0 )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTBMV ', 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 := A*X. * IF( LSAME( UPLO, 'U' ) )THEN KPLUS1 = K + 1 IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = X( J ) L = KPLUS1 - J DO 10, I = MAX( 1, J - K ), J - 1 X( I ) = X( I ) + TEMP*A( L + I, J ) 10 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( KPLUS1, J ) END IF 20 CONTINUE ELSE JX = KX DO 40, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX L = KPLUS1 - J DO 30, I = MAX( 1, J - K ), J - 1 X( IX ) = X( IX ) + TEMP*A( L + I, J ) IX = IX + INCX 30 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( KPLUS1, J ) END IF JX = JX + INCX IF( J.GT.K ) $ KX = KX + INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN TEMP = X( J ) L = 1 - J DO 50, I = MIN( N, J + K ), J + 1, -1 X( I ) = X( I ) + TEMP*A( L + I, J ) 50 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*A( 1, J ) END IF 60 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX L = 1 - J DO 70, I = MIN( N, J + K ), J + 1, -1 X( IX ) = X( IX ) + TEMP*A( L + I, J ) IX = IX - INCX 70 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*A( 1, J ) END IF JX = JX - INCX IF( ( N - J ).GE.K ) $ KX = KX - INCX 80 CONTINUE END IF END IF ELSE * * FORM X := A'*X. * IF( LSAME( UPLO, 'U' ) )THEN KPLUS1 = K + 1 IF( INCX.EQ.1 )THEN DO 100, J = N, 1, -1 TEMP = X( J ) L = KPLUS1 - J IF( NOUNIT ) $ TEMP = TEMP*A( KPLUS1, J ) DO 90, I = J - 1, MAX( 1, J - K ), -1 TEMP = TEMP + A( L + I, J )*X( I ) 90 CONTINUE X( J ) = TEMP 100 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 120, J = N, 1, -1 TEMP = X( JX ) KX = KX - INCX IX = KX L = KPLUS1 - J IF( NOUNIT ) $ TEMP = TEMP*A( KPLUS1, J ) DO 110, I = J - 1, MAX( 1, J - K ), -1 TEMP = TEMP + A( L + I, J )*X( IX ) IX = IX - INCX 110 CONTINUE X( JX ) = TEMP JX = JX - INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = 1, N TEMP = X( J ) L = 1 - J IF( NOUNIT ) $ TEMP = TEMP*A( 1, J ) DO 130, I = J + 1, MIN( N, J + K ) TEMP = TEMP + A( L + I, J )*X( I ) 130 CONTINUE X( J ) = TEMP 140 CONTINUE ELSE JX = KX DO 160, J = 1, N TEMP = X( JX ) KX = KX + INCX IX = KX L = 1 - J IF( NOUNIT ) $ TEMP = TEMP*A( 1, J ) DO 150, I = J + 1, MIN( N, J + K ) TEMP = TEMP + A( L + I, J )*X( IX ) IX = IX + INCX 150 CONTINUE X( JX ) = TEMP JX = JX + INCX 160 CONTINUE END IF END IF END IF * RETURN * * END OF DTBMV . * END 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 .