C************************************************************************* C * C * C SPARSE SVD VIA EIGENSYSTEM OF EQUIVALENT SHIFTED 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 TMS1 C.... C.... TMS1 IMPLEMENTS A TRACE-MINIMIZATION SVD METHOD FOR DETERMINING C.... THE P-LARGEST SINGULAR TRIPLETS OF A LARGE SPARSE MATRIX A. THIS C.... IS A PARALLEL METHOD WHICH PERMITS CONCURRENT ITERATIONS OF THE C.... CLASSICAL CONJUGATE GRADIENT METHOD. PARALLELIZATION DIRECTIVES C.... SUCH AS THE C$DOACROSS FOR THE SEQUENT SYMMETRY MULTIPROCESSOR C.... (SHARED MEMORY) PRECEDE LOOPS WHOSE ITERATIONS CAN BE INDEPENDENTLY C.... SCHEDULED ACROSS PROCESSORS OF ANY PARALLEL MACHINE. THESE C.... DIRECTIVES WILL BE TREATED AS COMMENTS BY ANY OTHER MACHINE'S C.... PREPROCESSOR, OF COURSE, AND ITERATIONS OF THE CORRESPONDING LOOPS C.... WILL BE EXECUTED IN SEQUENTIAL FASHION. C.... C.... SUBROUTINE TSVD1 COMPUTES THE SINGULAR TRIPLETS OF THE MATRIX A C.... VIA THE EQUIVALENT SYMMETRIC EIGENSYSTEM OF C.... C.... [ALPHA*I A] C.... B = [ ] , WHERE A IS M (=NROW) BY N (=NCOL), C.... [A' ALPHA*I] C.... C.... SO THAT {U,ABS(LAMBDA-ALPHA),V} IS A SINGULAR TRIPLET OF A. C.... ALPHA IS CHOSEN SO THAT B IS (SYMMETRIC) POSITIVE DEFINITE. IN C.... THE CALLING PROGRAM, ALPHA IS DETERMINED AS THE MATRIX 1-NORM OF C.... A. THE USER CAN SET ALPHA TO BE ANY KNOWN UPPER BOUND FOR THE C.... THE LARGEST SINGULAR VALUES OF THE MATRIX A. C.... (A' = TRANSPOSE OF A) C C.... USER SUPPLIED ROUTINES: OPAT, OPB C.... C.... OPAT(X,Y) TAKES AN M-VECTOR X AND SHOULD RETURN A'*X IN Y. C.... OPB(N,X,Y,SHIFT) TAKES AN N-VECTOR X AND SHOULD RETURN D*X C.... IN Y, WHERE D=[B-SHIFT*I], I IS THE N-TH ORDER IDENTITY MATRIX. C.... C.... USER SHOULD REPLACE CALLS TO DTIME WITH AN APPROPRIATE CALL TO C.... AN INSTRINSIC TIMING ROUTINE THAT RETURNS DELTA USER CPU TIME C.... (I.E., USER CPU TIME ELAPSED FROM PREVIOUS CALL TO TIMING ROUTINE). C.... ELAPSE RETURNS ELAPSED WALL-CLOCK TIME (MACHINE DEPENDENT) C.... C.... TMS1 UTILIZES RITZ SHIFTING FOR ACCELERATION OF CONVERGENCE. C.... C.....PARAMETERS IN TMS1 DRIVER: C.... C.... LDY SHOULD BE AT LEAST AS LARGE AS THE ORDER OF MATRIX B C.... ISMAX IS THE MAXIMUM INITIAL SUBSPACE DIMENSION FOR SUBR. TSVD1 C.... MAXI IS THE MAXIMUM NUMBER OF TRACE MIN. ITERATIONS ALLOWED. C.... C.....PARAMETERS IN TMS1 DRIVER AND SUBROUTINE TSVD1: C.... C.... NMAX IS THE MAXIMUM COLUMN DIMENSION FOR THE SPARSE MATRIX A C.... NZMAX IS THE MAXIMUM NUMBER OF NONZEROS FOR THE SPARSE MATRIX A C.... C.... PARAMETER IN SUBROUTINES CGT AND CGTS: C.... C.... NCG SHOULD BE AT LEAST AS LARGE AS THE ORDER OF MATRIX B. C.... PARAMETER(LDY=500,ISMAX=12,LDW1=LDY,LDW2=LDW1,LDW3=ISMAX, * LDW4=LDW1,LDW5=ISMAX,LDWI=ISMAX, * MAXI=150,NMAX=500,NZMAX=2000) C DOUBLE PRECISION WORK1(LDW1,ISMAX), WORK2(LDW2,ISMAX), * WORK3(LDW3,ISMAX), WORK4(LDW4,3), * WORK5(LDW5,4), Y(LDY,MAXI), * RES(MAXI),SIG(MAXI), * VALUE(NZMAX), ALPHA, TOL, RED, DSUM C REAL TIME(5),WCLK0,WCLK1 C INTEGER N, P, S, JOB, I, II, J, JJ, LDY, * TITER(MAXI), MEM, ISMAX, NMAX, NZMAX, * ITCG(MAXI), MAXI, MXV(3), ITIME(4), * IWORK(LDWI,2), POINTR(NMAX), ROWIND(NZMAX), * NCOL, NROW, MATRIX, RESULT, INPUT, * NNZERO, NRHS, IVEC, IWALL, NINT, * LDWI, LDW1, LDW2, LDW3, LDW4, LDW5, * IERR LOGICAL LWORK(ISMAX),VEC C C------------------------------------------------------------ C C ----- HARWELL-BOEING COLLECTION DATA STRUCTURE C C------------------------------------------------------------ C CHARACTER TITLE*72, KEY*8, TYPE*3, PTRFMT*16, INDFMT*16, * VALFMT*20, RHSFMT*20 C COMMON /SPARSE/ VALUE,ALPHA,POINTR,ROWIND,NCOL,NROW C.... INPUT=5 MATRIX=7 RESULT=9 C.... C.... OPEN FILES FOR INPUT/OUTPUT C.... OPEN(MATRIX,FILE='LAI7') OPEN(RESULT,FILE='TMO9') OPEN(INPUT,FILE='TMI5') C.... C.... READ IN PARAMETERS FOR THIS TEST RUN C.... C.... P := NUMBER OF SINGULAR TRIPLETS (LARGEST) DESIRED C.... S := INITIAL SUBSPACE DIMENSION C.... JOB := CONTROLS USE OF RITZ SHIFTING (0=NO,1=YES) C.... TOL := USER-SPECIFIED TOLERANCE FOR RESIDUALS C.... RED := RESIDUAL NORM REDUCTION FACTOR FOR INITIATING SHIFTING C.... VEC := OUTPUT SINGULAR VECTORS? (BOOLEAN) C.... READ(INPUT,*) P, S, JOB, TOL, RED, VEC C.... IF (VEC) THEN IVEC=8 OPEN(IVEC,FILE='TMO8',FORM='UNFORMATTED') ENDIF 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 TOO BIG ---- C IF (NCOL .GE. NMAX .OR. NNZERO .GT. NZMAX) THEN WRITE(RESULT,*) ' 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 N=NROW+NCOL C C.... ESTIMATE ALPHA (VIA MATRIX 1-NORM) C.... (USED FOR UPPER BOUND ON MAX SINGULAR VALUE OF MATRIX A) C ALPHA=0.0D0 DO 75 JJ=1,NCOL DSUM=0.0D0 DO 85 II=POINTR(JJ),POINTR(JJ+1)-1 DSUM=DSUM+DABS(VALUE(II)) 85 CONTINUE ALPHA=DMAX1(ALPHA,DSUM) 75 CONTINUE C CALL ELAPSE(WCLK0) CALL TSVD1 (N,P,S,JOB, 1 MAXI,TOL,RED,SIG,Y,LDY,MEM, 2 TITER,ITCG,TIME,RES,MXV,WORK1,LDW1, 3 WORK2,LDW2,WORK3,LDW3,WORK4,LDW4, 4 WORK5,LDW5,IWORK,LDWI,LWORK,IERR) CALL ELAPSE(WCLK1) IWALL=NINT(MOD(WCLK1-WCLK0+1000000.0,1000000.0)) ITIME(1)=NINT(100*(TIME(1)/TIME(5))) ITIME(2)=NINT(100*(TIME(2)/TIME(5))) ITIME(3)=NINT(100*(TIME(3)/TIME(5))) ITIME(4)=NINT(100*(TIME(4)/TIME(5))) C WRITE(RESULT,2000) N,MAXI,TITER(P),P,S,MEM,JOB,VEC,ALPHA, * TOL,RED,MXV(1),MXV(2),MXV(3),IERR WRITE(RESULT,2001) TIME(1),ITIME(1), * TIME(2),ITIME(2), * TIME(3),ITIME(3), * TIME(4),ITIME(4), * TIME(5),IWALL,TITLE,KEY,NROW,NCOL,N 2000 FORMAT( * 1X,'... ' * /1X,'... TRACE MINIMIZATION ON THE ' * /1X,'... EQUIVALENT CYCLIC E-PROBLEM' * /1X,'... NO. OF EQUATIONS =',I10 * /1X,'... MAX. NO. OF ITERATIONS =',I10 * /1X,'... NO. OF ITERATIONS TAKEN =',I10 * /1X,'... NO. OF DESIRED EIGENPAIRS =',I10 * /1X,'... INITIAL SUBSPACE DIM. =',I10 * /1X,'... MEMORY REQUIRED (BYTES) =',I10 * /1X,'... JOB: 0=NO SHIFT, 1=SHIFT =',I10 * /1X,'... WANT S-VECTORS? [T/F] =',L10 * /1X,'... ALPHA =',1PE10.2 * /1X,'... FINAL RESIDUAL TOLERANCE =',1PE10.2 * /1X,'... RESIDUAL REDUCTION TOL. =',1PE10.2 * /1X,'... MULTIPLICATIONS BY A =',I10 * /1X,'... MULTIPLICATIONS BY A^T =',I10 * /1X,'... TOTAL NUMBER OF MULT. =',I10 * /1X,'... IERR FROM SUBR. TSVD1 =',I10 * /1X,'... ') 2001 FORMAT( * 1X,'... CPU TIME BREAKDOWN: ', * /1X,'... GRAM-SCHMIDT ORTHOG. =',1PE10.2,2X,'(',I3,'%)' * /1X,'... SPECTRAL DECOMPOSITION =',1PE10.2,2X,'(',I3,'%)' * /1X,'... CONVERGENCE CRITERIA =',1PE10.2,2X,'(',I3,'%)' * /1X,'... CONJUGATE GRADIENT =',1PE10.2,2X,'(',I3,'%)' * /1X,'... TOTAL CPU TIME (SEC) =',1PE10.2 * /1X,'... WALL-CLOCK TIME (SEC) =',I10 * //1X,A50 * / 1X,A50 * /1X,'... NO. OF DOCUMENTS (ROWS) = ',I8 * /1X,'... NO. OF TERMS (COLS) = ',I8 * /1X,'... ORDER OF MATRIX B = ',I8 * /1X,'... ') C WRITE(RESULT,9999) (I,SIG(I),RES(I),TITER(I),ITCG(I), * I=1,P) IF (VEC) THEN DO 35 J=1,P WRITE(IVEC) (Y(I,J),I=1,N) 35 CONTINUE ENDIF C 9999 FORMAT(1X,'...... ' * /1X,'...... ',4X,'COMPUTED SINGULAR VALUES',2X, * '(','RESIDUAL NORMS',')',2X,'T-MIN STEPS',3X, * 'CG STEPS', * /1X,'...... ' * /(1X,'...... ',I3,3X,1PE22.14,2X,'(',1PE11.2,')',2I12)) STOP END C C ELAPSE RETURNS WALL CLOCK TIME, T, IN SECONDS. C NOTE T IS OF TYPE 'REAL' C C THE INTRINSIC FUNTION 'TIME' IS CALLED. C TIME RETURNS SYSTEM TIME IN SECONDS SINCE JAN. 1, 1970, C AND WE ONLY USE THE LEAST SIGNIFICANT 6 DIGITS. C SUBROUTINE ELAPSE(T) REAL T INTEGER TIME T=MOD(TIME(),1000000) RETURN END C SUBROUTINE TSVD1(N,P,S,JOB,MAXI,TOL,RED,SIG,Y,LDY,MEM, * TITER,ITCGT,TIME,RES,MXV,WORK1,LDW1, * WORK2,LDW2,WORK3,LDW3,WORK4,LDW4,WORK5, * LDW5,IWORK,LDWI,LWORK,IERR) C C.... TSVD1 IMPLEMENTS A TRACE-MINIMIZATION SVD METHOD FOR DETERMINING C.... THE P-LARGEST SINGULAR TRIPLETS OF A LARGE SPARSE MATRIX A. THIS C.... IS A PARALLEL METHOD WHICH PERMITS CONCURRENT ITERATIONS OF THE C.... CLASSICAL CONJUGATE GRADIENT METHOD. PARALLELIZATION DIRECTIVES C.... SUCH AS THE C$DOACROSS FOR THE SEQUENT SYMMETRY MULTIPROCESSOR C.... (SHARED MEMORY) PRECEDE LOOPS WHOSE ITERATIONS CAN BE INDEPENDENTLY C.... SCHEDULED ACROSS PROCESSORS OF ANY PARALLEL MACHINE. THESE C.... DIRECTIVES WILL BE TREATED AS COMMENTS BY ANY OTHER MACHINE'S C.... PREPROCESSOR, OF COURSE, AND ITERATIONS OF THE CORRESPONDING LOOPS C.... WILL BE EXECUTED IN SEQUENTIAL FASHION. C.... C.... TMS1 COMPUTES THE SINGULAR TRIPLETS OF THE MATRIX A VIA THE C.... EQUIVALENT SYMMETRIC EIGENVALUE PROBLEM C.... C.... [ALPHA*I A] C.... B = [ ] , WHERE A IS M (=NROW) BY N (=NCOL), C.... [A' ALPHA*I] C.... C.... SO THAT {U,ABS(LAMBDA-ALPHA),V} IS A SINGULAR TRIPLET OF A. C.... ALPHA IS CHOSEN SO THAT B IS (SYMMETRIC) POSITIVE DEFINITE. IN C.... THE CALLING PROGRAM, ALPHA IS DETERMINED AS THE MATRIX 1-NORM OF C.... A. THE USER CAN SET ALPHA TO BE ANY KNOWN UPPER BOUND FOR THE C.... THE LARGEST SINGULAR VALUES OF THE MATRIX A. C.... (A' = TRANSPOSE OF A) C C.... USER SUPPLIED ROUTINES: OPAT, OPB C.... C.... OPAT(X,Y) TAKES AN M-VECTOR X AND SHOULD RETURN A'*X IN Y. C.... OPB(NROW+NCOL,X,Y,SHIFT) TAKES AN (M+N)-VECTOR X AND SHOULD RETURN D*X C.... IN Y, WHERE D=[B-SHIFT*I], I IS THE (M+N) ORDER IDENTITY MATRIX. C.... C.... USER SHOULD REPLACE CALLS TO DTIME WITH AN APPROPRIATE CALL TO C.... AN INSTRINSIC TIMING ROUTINE THAT RETURNS DELTA USER CPU TIME C.... (I.E., USER CPU TIME ELAPSED FROM PREVIOUS CALL TO TIMING ROUTINE). C.... C.... TSVD1 UTILIZES RITZ SHIFTING FOR ACCELERATION OF CONVERGENCE. C.... C.... INPUT PARAMETERS: C.... C.... N (INTEGER) ORDER OF EQUIVALENT SYMMETRIC EIGENVALUE PROBLEM. C.... SHOULD BE EQUAL TO NROWS+NCOLS, WHERE NROWS IS THE C.... NUMBER OF ROWS OF A AND NCOLS IS THE NUMBER OF C.... COLUMNS OF A. C.... C.... P (INTEGER) NUMBER OF DESIRED SINGULAR TRIPLETS (LARGEST) OF A. C.... C.... S (INTEGER) DIMENSION OF INITIAL SUBSPACE (S SHOULD BE GREATER C.... THAN P; S=2*P IS USUALLY SAFE BUT SINCE THE COMPLEXITY C.... OF THE METHOD IS DETERMINED BY S, THE USER SHOULD C.... TRY TO KEEP S AS CLOSE TO P AS POSSIBLE). C.... C.... JOB (INTEGER) ACCELERATION STRATEGY SWITCH: C.... C.... JOB = 0, NO ACCELERATION USED, C.... = 1, RITZ-SHIFTING USED. C.... C.... MAXI (INTEGER) MAXIMUM NUMBER OF TRACE MINIMIZATION STEPS ALLOWED. C.... C.... TOL (DOUBLE USER-SUPPLIED TOLERANCE FOR RESIDUAL OF AN C.... PRECISION) EQUIVALENT EIGENPAIR OF B (SINGULAR TRIPLET OF A). C.... C.... RED (DOUBLE USER-SUPPLIED TOLERANCE FOR RESIDUAL REDUCTION TO C.... PRECISION) INITIATE RITZ-SHIFTING WHEN JOB=1. C.... C.... LDY (INTEGER) LEADING DIMENSION OF ARRAY Y, IT SHOULD BE GREATER C.... THAN OR EQUAL TO NROWS+NCOLS (ORDER OF MATRIX B). C.... C.... LDW1 (INTEGER) LEADING DIMENSION OF ARRAY WORK1, IT SHOULD BE GREATER C.... THAN OR EQUAL TO LDY. C.... C.... WORK1(LDW1,S) (DOUBLE PRECISION) WORK ARRAY. C.... C.... LDW2 (INTEGER) LEADING DIMENSION OF ARRAY WORK2, IT SHOULD BE GREATER C.... THAN OR EQUAL TO LDY. C.... C.... WORK2(LDW2,S) (DOUBLE PRECISION) WORK ARRAY. C.... C.... LDW3 (INTEGER) LEADING DIMENSION OF ARRAY WORK3, IT SHOULD BE GREATER C.... THAN OR EQUAL TO S. C.... C.... WORK3(LDW3,S) (DOUBLE PRECISION) WORK ARRAY. C.... C.... LDW4 (INTEGER) LEADING DIMENSION OF ARRAY WORK4, IT SHOULD BE GREATER C.... THAN OR EQUAL TO LDY. C.... C.... WORK4(LDW4,3) (DOUBLE PRECISION) WORK ARRAY. C.... C.... LDW5 (INTEGER) LEADING DIMENSION OF ARRAY WORK5, IT SHOULD BE GREATER C.... THAN OR EQUAL TO S. C.... C.... WORK5(LDW5,4) (DOUBLE PRECISION) WORK ARRAY. C.... C.... LDWI (INTEGER) LEADING DIMENSION OF ARRAY IWORK, IT SHOULD BE GREATER C.... THAN OR EQUAL TO S. C.... C.... IWORK(LDWI,2) (INTEGER) WORK ARRAY. C.... C.... LWORK(S) (LOGICAL) WORK ARRAY. C.... C.... OUTPUT PARAMETERS: C.... C.... IERR (INTEGER) ERROR FLAG: C.... IERR=99, INPUT PARAMETER JOB INVALID. C.... C.... MEM (INTEGER) NUMBER OF BYTES OF MEMORY NEEDED C.... (REAL, INTEGER, LOGICAL = 4 BYTES, C.... DOUBLE PRECISION = 8 BYTES) C.... C.... MXV(3) (INTEGER) MATRIX TIMES VECTOR COUNTERS: C.... MXV(1) = NUMBER OF A*X, C.... MXV(2) = NUMBER OF A'*X, C.... MXV(3) = MXV(1)+MXV(2) C.... C.... SIG(P) (DOUBLE PRECISION) CONTAINS P-DESIRED SINGULAR VALUES. C.... C.... Y(LDY,P) (DOUBLE PRECISION) CONTAINS LEFT AND RIGHT SINGULAR C.... VECTORS, I.E., C.... C.... Y(1:N,1:P) = [U(1:R,1:P)' | V(1:C,1:P)']', C.... WHERE R = NO. OF ROWS OF MATRIX A AND C.... C = NO. OF COLS OF MATRIX A. C.... C.... C.... TITER(P) (INTEGER) TITER(I) := NUMBER OF TRACE MIN. STEPS C.... NEED FOR I-TH TRIPLET OF A. C.... C.... C.... ITCGT(S) (INTEGER) ITCGT(I) := NUMBER OF CG STEPS NEEDED FOR C.... I-TH TRIPLET OF A. (ARRAY MUST HAVE S ELEMENTS) C.... C.... TIME(5) (REAL) TIMING BREAKDOWN ARRAY: C.... C.... TIME(1) = GRAM-SCHMIDT ORTHOGONALIZATION C.... TIME(2) = SECTION-FORMATION (SPECTRAL DECOMPOSITION) C.... TIME(3) = CONVERGENCE CRITERIA C.... TIME(4) = CONJUGATE GRADIENT METHOD C.... TIME(5) = TOTAL EXECUTION TIME (USER-CPU) C.... C.... RES(P) (DOUBLE PRECISION) 2-NORMS OF RESIDUAL VECTORS C.... (A*Y[I]-SIG(I)*Y[I]), I=1,2,...,P. C.... C.... C.... ROUTINES CALLED: C.... C.... (TIMER) DTIME C.... (RNG) RANDOM C.... (PRECISION) EPSLON C.... (MATH) DSQRT,DABS,DMAX1,MIN0 C.... (BLAS) DAXPY,DDOT,DNRM2 C.... (BLAS3) DGEMM,DGEMV,XERBLA C.... (EISPACK) TRED2,TQL2,PYTHAG C.... C.... INTEGER N, P, LDY, S, MEM, JOB, IPTR, LEFT, * I, J, II, JJ, K,IERR, IRAND, ITER, * MAXI, TITER(P), ITCGT(S), MXV(3), * LDW1, LDW2, LDW3, LDW4, LDW5, LDWI, * IWORK(LDWI,2), NMAX, NZMAX REAL T1, TOTAL, DT(2), DTIME, * SEC1, SEC2, SEC21, SEC22, SEC23, * SEC3, SEC4, TIME(5) DOUBLE PRECISION WORK1(LDW1,S), WORK2(LDW2,S), WORK3(LDW3,S), * WORK4(LDW4,3), WORK5(LDW5,4), * Y(LDY,S), SIG(S), RES(P), EPSLON, MEPS, * DDOT, DNRM2, DASUM, RANDOM, RED, TOL LOGICAL LWORK(S),SHIFT C C.... COMMON BLOCKS: C PARAMETER(NMAX=500,NZMAX=2000) DOUBLE PRECISION VALUE(NZMAX), ALPHA INTEGER POINTR(NMAX), ROWIND(NZMAX), NCOL, NROW COMMON /SPARSE/ VALUE,ALPHA,POINTR,ROWIND,NCOL,NROW C C.....GET BYTE COUNT (ESTIMATE) C MEM = 4*(48+2*P+2*LDWI)+ * 8*(4+P+S+3*LDW4+4*LDW5+ * S*(LDW1+LDW2+LDW3+LDY) ) C C.... GET MACHINE EPSILON (MEPS) C MEPS=EPSLON(1.0D0) SHIFT =.FALSE. IERR=0 IF(JOB.NE.0.AND.JOB.NE.1) THEN IERR=99 RETURN ENDIF C IF(JOB.EQ.1) THEN DO 36 I=1,P LWORK(I) =.FALSE. WORK5(I,3)=-1.0D0 WORK5(I,4)=-1.0D0 36 CONTINUE ENDIF C C.... INITIALIZE TIMERS AND COUNTERS: C SEC1 = 0.0 SEC21 = 0.0 SEC22 = 0.0 SEC23 = 0.0 SEC3 = 0.0 SEC4 = 0.0 MXV(1)= 0 MXV(2)= 0 C C-- INITIALIZE Y(1:N,1:S) = RANDOM MATRIX C C-- (CARRY S VECTORS IN THE TMIN ITERATIONS, ASSUMING S.GE.P) C IRAND = 91827211 C DO 30 K = 1, S SIG(K) = 0.0D0 WORK5(K,2) = 0.0D0 30 CONTINUE C DO 31 K = 1, S DO 40 I = 1, N Y(I,K) = RANDOM(IRAND) 40 CONTINUE 31 CONTINUE C C-------------------------------------------------------------- C C POINTER AND COUNTER FOR HYBRID MONITOR: C C 1 2 3 4 5 ... I ... P C ----------------------- C SIG:| | | | | |...| |...| | (ASCENDING ORDER) C ----------------------- C ^ C | C IPTR : POINTS TO FIRST E-VALUE OF B C THAT HAS NOT CONVERGED C C LEFT:=S-IPTR+1 (IE HOW MANY Y-VECTORS REMAIN FOR TSVD1) C-------------------------------------------------------------- C C INITIALIZE A FEW POINTERS AND COUNTERS: C TOTAL=0.0 IPTR=1 LEFT=S DO 170 II=1,P TITER(II)=0 170 CONTINUE DO 175 II=1,S ITCGT(II)=0 175 CONTINUE C C-------------------------------------------------------------- C MAIN TMIN ITERATION LOOP (NMAX ITERATIONS) C-------------------------------------------------------------- C DO 500 ITER = 1,MAXI C C.... SECTION FORMATION C C C.... USE MODIFIED GRAM-SCHMIDT FOR B-ORTHOGONALIZATION C T1 = DTIME(DT(1)) DO 615 JJ=1,S DO 616 II=1,N WORK2(II,JJ)=Y(II,JJ) 616 CONTINUE 615 CONTINUE CALL ORTHG(N,0,S,WORK1,LDW1,WORK2,LDW2,WORK4(1,1)) T1 = DTIME(DT(1)) C SEC1 = SEC1 + T1 C C C.... FORM WORK1(1:N,IPTR:S) = B(1:N,1:N)*WORK2(1:N,IPTR:S) C T1 = DTIME(DT(1)) C C.... ECONOMIZATION OF 2-CYCLIC FORM USED HERE C DO 240 II=IPTR,S CALL OPAT(WORK2(1,II),WORK1(1,II)) 240 CONTINUE MXV(2)=MXV(2)+S-IPTR+1 C C.... FORM WORK3(1:LEFT,1:LEFT) = WORK2(1:N,IPTR:S)'*WORK1(1:N,IPTR:S) C CALL DGEMM('T','N',LEFT,LEFT,NCOL,1.0D0,WORK2(NROW+1, * IPTR),LDW2,WORK1(1,IPTR),LDW1,0.0D0,WORK3,LDW3) C DO 1620 JJ=1,LEFT DO 1610 II=JJ,LEFT WORK3(II,JJ)=WORK3(II,JJ)+WORK3(JJ,II) 1610 CONTINUE DO 1615 II=JJ,LEFT WORK3(JJ,II)=WORK3(II,JJ) 1615 CONTINUE WORK3(JJ,JJ)=WORK3(JJ,JJ)+ALPHA 1620 CONTINUE C C.... LOAD GERSHGORIN RADII C IF(JOB.EQ.1) THEN DO 1630 J = 1,LEFT WORK4(J,1)=DASUM(LEFT,WORK3(J,1),LDW3)-DABS(WORK3(J,J)) 1630 CONTINUE ENDIF C C.... COMPUTE THE EIGENVALUES AND EIGENVECTORS OF WORK3: C T1 = DTIME(DT(1)) SEC21 = SEC21 + T1 C T1 = DTIME(DT(1)) C C.... EIGENVECTORS OVERWRITE ARRAY WORK3(:,:) C.... STORE CURRENT E-VALUES IN WORK5(:,2) C DO 560 JJ=IPTR,S WORK5(JJ,2)=SIG(JJ) 560 CONTINUE C CALL TRED2(LDW3,LEFT,WORK3,SIG(IPTR),WORK4(1,2),WORK3) CALL TQL2(LDW3,LEFT,SIG(IPTR),WORK4(1,2),WORK3,IERR) C T1=DTIME(DT(1)) C SEC22 = SEC22 + T1 C C.... FORM Y(1:N,IPTR:S) = WORK2(1:N,IPTR:S)*WORK3(1:LEFT,1:LEFT) C T1=DTIME(DT(1)) C CALL DGEMM('N','N',N,LEFT,LEFT,1.0D0,WORK2(1,IPTR), * LDW2,WORK3,LDW3,0.0D0,Y(1,IPTR),LDY) C T1=DTIME(DT(1)) SEC23 = SEC23 + T1 C C.... TEST FOR CONVERGENCE HERE C T1=DTIME(DT(1)) C DO 840 J=IPTR,S CALL OPB(N,Y(1,J),WORK2(1,J),0.0D0) CALL DAXPY(N,-SIG(J),Y(1,J),1,WORK2(1,J),1) WORK4(J-IPTR+1,3) = DNRM2(N,WORK2(1,J),1) IF(J.LE.P) TITER(J)=TITER(J)+1 840 CONTINUE MXV(1) =MXV(1)+S-IPTR+1 MXV(2) =MXV(2)+S-IPTR+1 C T1=DTIME(DT(1)) SEC3 = SEC3 + T1 C C.... ARRAY WORK4(:,3) STORES THE VECTOR 2-NORMS OF RESIDUAL VECTORS C DO 50 K=1,LEFT IF (DABS(WORK4(K,3)).LE.TOL) THEN IPTR=IPTR+1 IF (IPTR.GT.P) GO TO 400 LEFT=S-IPTR+1 ENDIF 50 CONTINUE C IF(.NOT.SHIFT.AND.JOB.EQ.1) THEN C C.... WORK4(:,1) STORES GERSHGORIN RADII C CALL DISK(LEFT,SIG(IPTR),WORK4(1,1),IWORK(1,1), * IWORK(1,2)) C C.... IWORK(:,1) IS THE CLUSTERING VECTOR C DO 244 K = IPTR,S C IF(IWORK(K-IPTR+1,1).EQ.1) THEN CALL ISOL(K,WORK4(K-IPTR+1,3), * RED,LWORK,WORK5(1,3),WORK5(1,4),S) ELSE IF(IWORK(K-IPTR+1,1).GT.1) THEN CALL CLUS(K,IWORK(K-IPTR+1,1), * WORK4(K-IPTR+1,3),RED, * LWORK,WORK5(1,3),WORK5(1,4),WORK4(1,2),S) ENDIF C 244 CONTINUE C ENDIF C C CONTINUE ALGORITHM... C C.... SHIFT START C IF(ITER.GT.1.AND..NOT.SHIFT.AND.JOB.EQ.1) THEN C IF(LWORK(IPTR)) THEN SHIFT=.TRUE. CALL ORTHG(N,0,S,WORK1,LDW1,Y,LDY,WORK4(1,1)) ENDIF ENDIF C IF(SHIFT) THEN C C.... COMPUTE SHIFTS IN WORK5(:,1) HERE: C DO 638 K=IPTR,S IF(K.GT.IPTR.AND.K.LE.P) THEN WORK5(K,1)=0.0 C DO 639 J=IPTR-1,K-1 IF(J.NE.0.AND.SIG(J).LE.SIG(K) * -WORK4(K-IPTR+1,3)) * WORK5(K,1)=SIG(J) 639 CONTINUE C IF(WORK5(K,1).EQ.SIG(K-1).AND. * WORK5(K-1,1).EQ.SIG(K-1)) * WORK5(K,1)=SIG(K) IF(WORK5(K,1).EQ.0.0) * WORK5(K,1)=WORK5(K-1,1) ELSE IF(K.GT.P) THEN WORK5(K,1)=WORK5(P,1) ELSE IF(K.EQ.IPTR) THEN WORK5(K,1)=SIG(K) ENDIF 638 CONTINUE ENDIF C T1 = DTIME(DT(1)) C C.... DO PARALLEL CG ITERATIONS C C.... LOAD Y INTO WORK1 ARRAY FOR ORTHOG. PROJECTOR IN SUBROUTINE CGT C DO 490 JJ=1,LEFT DO 495 II=1,N WORK1(II,JJ)=Y(II,IPTR+JJ-1) 495 CONTINUE 490 CONTINUE C C.... CG LOOP FOR INDEPENDENT SYSTEMS (NO SHIFTING USED) C IF(.NOT.SHIFT) THEN C$DOACROSS SHARE(Y,LDW1,ITCGT,SIG,IWORK,MEPS,N,WORK1, C$ * IPTR,S,LEFT) C$ * LOCAL(II) DO 515 II=1,LEFT CALL CGT (N,LEFT,Y(1,IPTR+II-1),WORK1,LDW1, * ITCGT(IPTR+II-1),SIG(IPTR+II-1),SIG(S), * IWORK(II,1),MEPS) C 515 CONTINUE C ELSE C C.... SHIFT WITH WORK5(:,1) IN CG ITERATIONS C C$DOACROSS SHARE(Y,LDW1,ITCGT,SIG,WORK5,IWORK,MEPS,N, C$ * WORK1,IPTR,S,LEFT) C$ * LOCAL(II) DO 516 II=1,LEFT CALL CGTS (N,LEFT,Y(1,IPTR+II-1),WORK1,LDW1, * ITCGT(IPTR+II-1),SIG(IPTR+II-1), * WORK5(IPTR+II-1,2),SIG(S),WORK5(IPTR+II-1,1), * IWORK(II,1),MEPS) 516 CONTINUE C ENDIF DO 517 II=1,LEFT MXV(1)=MXV(1)+IWORK(II,1) MXV(2)=MXV(2)+IWORK(II,1) 517 CONTINUE C T1 = DTIME(DT(1)) C SEC4 = SEC4 + T1 C 500 CONTINUE 400 CONTINUE C C.... COMPUTE FINAL 2-NORMS OF RESIDUAL VECTORS W.R.T. B (INDEFINITE) C DO 740 J = 1,P CALL OPB(N,Y(1,J),WORK2(1,J),ALPHA) SIG(J)=DDOT(N,Y(1,J),1,WORK2(1,J),1) * /DDOT(N,Y(1,J),1, Y(1,J),1) CALL DAXPY(N,-SIG(J),Y(1,J),1,WORK2(1,J),1) WORK4(J,3) = DNRM2(N,WORK2(1,J),1) SIG(J)=DABS(SIG(J)) 740 CONTINUE C MXV(1) = MXV(1)+P MXV(2) = MXV(2)+P C SEC2 = SEC21 + SEC22 + SEC23 TOTAL = TOTAL + SEC1 + SEC2 + SEC3 + SEC4 C C.... LOAD OUTPUT TIME AND MXV ARRAYS C TIME(1)=SEC1 TIME(2)=SEC2 TIME(3)=SEC3 TIME(4)=SEC4 TIME(5)=TOTAL MXV(3)=MXV(1)+MXV(2) C C.... LOAD RESIDUAL VECTOR C DO 430 II=1,P RES(II)=WORK4(II,3) 430 CONTINUE C RETURN END C C.... END OF TSVD1 C DOUBLE PRECISION FUNCTION EPSLON (X) C C ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. C DOUBLE PRECISION A, B, C, EPS, X C C THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS C SATISFYING THE FOLLOWING TWO ASSUMPTIONS, C 1. THE BASE USED IN REPRESENTING FLOATING POINT C NUMBERS IS NOT A POWER OF THREE. C 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO C THE ACCURACY USED IN FLOATING POINT VARIABLES C THAT ARE STORED IN MEMORY. C THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO C FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING C ASSUMPTION 2. C UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT, C A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, C B HAS A ZERO FOR ITS LAST BIT OR DIGIT, C C IS NOT EXACTLY EQUAL TO ONE, C EPS MEASURES THE SEPARATION OF 1.0 FROM C THE NEXT LARGER FLOATING POINT NUMBER. C THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED C ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD. C C THIS VERSION DATED 4/6/83. C A = 4.0D0/3.0D0 10 B = A - 1.0D0 C = B + B + B EPS = DABS(C-1.0D0) IF (EPS .EQ. 0.0D0) GO TO 10 EPSLON = EPS*DABS(X) 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 DOUBLE PRECISION ALPHA, BETA * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * PURPOSE * ======= * * DGEMM PERFORMS ONE OF THE MATRIX-MATRIX OPERATIONS * * C := ALPHA*OP( A )*OP( B ) + BETA*C, * * WHERE OP( X ) IS ONE OF * * OP( X ) = X OR OP( X ) = X', * * ALPHA AND BETA ARE SCALARS, AND A, B AND C ARE MATRICES, WITH OP( A ) * AN M BY K MATRIX, OP( B ) A K BY N MATRIX AND C AN M BY N MATRIX. * * PARAMETERS * ========== * * TRANSA - CHARACTER*1. * ON ENTRY, TRANSA SPECIFIES THE FORM OF OP( A ) TO BE USED IN * THE MATRIX MULTIPLICATION AS FOLLOWS: * * TRANSA = 'N' OR 'N', OP( A ) = A. * * TRANSA = 'T' OR 'T', OP( A ) = A'. * * TRANSA = 'C' OR 'C', OP( A ) = A'. * * UNCHANGED ON EXIT. * * TRANSB - CHARACTER*1. * ON ENTRY, TRANSB SPECIFIES THE FORM OF OP( B ) TO BE USED IN * THE MATRIX MULTIPLICATION AS FOLLOWS: * * TRANSB = 'N' OR 'N', OP( B ) = B. * * TRANSB = 'T' OR 'T', OP( B ) = B'. * * TRANSB = 'C' OR 'C', OP( B ) = B'. * * UNCHANGED ON EXIT. * * M - INTEGER. * ON ENTRY, M SPECIFIES THE NUMBER OF ROWS OF THE MATRIX * OP( A ) AND OF THE MATRIX C. M MUST BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE NUMBER OF COLUMNS OF THE MATRIX * OP( B ) AND THE NUMBER OF COLUMNS OF THE MATRIX C. N MUST BE * AT LEAST ZERO. * UNCHANGED ON EXIT. * * K - INTEGER. * ON ENTRY, K SPECIFIES THE NUMBER OF COLUMNS OF THE MATRIX * OP( A ) AND THE NUMBER OF ROWS OF THE MATRIX OP( B ). K MUST * BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * ALPHA - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, KA ), WHERE KA IS * K WHEN TRANSA = 'N' OR 'N', AND IS M OTHERWISE. * BEFORE ENTRY WITH TRANSA = 'N' OR 'N', THE LEADING M BY K * PART OF THE ARRAY A MUST CONTAIN THE MATRIX A, OTHERWISE * THE LEADING K BY M PART OF THE ARRAY A MUST CONTAIN THE * MATRIX A. * UNCHANGED ON EXIT. * * LDA - INTEGER. * ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED * IN THE CALLING (SUB) PROGRAM. WHEN TRANSA = 'N' OR 'N' THEN * LDA MUST BE AT LEAST MAX( 1, M ), OTHERWISE LDA MUST BE AT * LEAST MAX( 1, K ). * UNCHANGED ON EXIT. * * B - DOUBLE PRECISION ARRAY OF DIMENSION ( LDB, KB ), WHERE KB IS * N WHEN TRANSB = 'N' OR 'N', AND IS K OTHERWISE. * BEFORE ENTRY WITH TRANSB = 'N' OR 'N', THE LEADING K BY N * PART OF THE ARRAY B MUST CONTAIN THE MATRIX B, OTHERWISE * THE LEADING N BY K PART OF THE ARRAY B MUST CONTAIN THE * MATRIX B. * UNCHANGED ON EXIT. * * LDB - INTEGER. * ON ENTRY, LDB SPECIFIES THE FIRST DIMENSION OF B AS DECLARED * IN THE CALLING (SUB) PROGRAM. WHEN TRANSB = 'N' OR 'N' THEN * LDB MUST BE AT LEAST MAX( 1, K ), OTHERWISE LDB MUST BE AT * LEAST MAX( 1, N ). * UNCHANGED ON EXIT. * * BETA - DOUBLE PRECISION. * ON ENTRY, BETA SPECIFIES THE SCALAR BETA. WHEN BETA IS * SUPPLIED AS ZERO THEN C NEED NOT BE SET ON INPUT. * UNCHANGED ON EXIT. * * C - DOUBLE PRECISION ARRAY OF DIMENSION ( LDC, N ). * BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY C MUST * CONTAIN THE MATRIX C, EXCEPT WHEN BETA IS ZERO, IN WHICH * CASE C NEED NOT BE SET ON ENTRY. * ON EXIT, THE ARRAY C IS OVERWRITTEN BY THE M BY N MATRIX * ( ALPHA*OP( A )*OP( B ) + BETA*C ). * * LDC - INTEGER. * ON ENTRY, LDC SPECIFIES THE FIRST DIMENSION OF C AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDC MUST BE AT LEAST * MAX( 1, M ). * UNCHANGED ON EXIT. * * * LEVEL 3 BLAS ROUTINE. * * -- WRITTEN ON 8-FEBRUARY-1989. * JACK DONGARRA, ARGONNE NATIONAL LABORATORY. * IAIN DUFF, AERE HARWELL. * JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD. * SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD. * * * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. LOCAL SCALARS .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB DOUBLE PRECISION TEMP * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. EXECUTABLE STATEMENTS .. * * SET NOTA AND NOTB AS TRUE IF A AND B RESPECTIVELY ARE NOT * TRANSPOSED AND SET NROWA, NCOLA AND NROWB AS THE NUMBER OF ROWS * AND COLUMNS OF A AND THE NUMBER OF ROWS OF B RESPECTIVELY. * NOTA = 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 C SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * DGEMV PERFORMS ONE OF THE MATRIX-VECTOR OPERATIONS * * Y := ALPHA*A*X + BETA*Y, OR Y := ALPHA*A'*X + BETA*Y, * * WHERE ALPHA AND BETA ARE SCALARS, X AND Y ARE VECTORS AND A IS AN * M BY N MATRIX. * * PARAMETERS * ========== * * TRANS - CHARACTER*1. * ON ENTRY, TRANS SPECIFIES THE OPERATION TO BE PERFORMED AS * FOLLOWS: * * TRANS = 'N' OR 'N' Y := ALPHA*A*X + BETA*Y. * * TRANS = 'T' OR 'T' Y := ALPHA*A'*X + BETA*Y. * * TRANS = 'C' OR 'C' Y := ALPHA*A'*X + BETA*Y. * * UNCHANGED ON EXIT. * * M - INTEGER. * ON ENTRY, M SPECIFIES THE NUMBER OF ROWS OF THE MATRIX A. * M MUST BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE NUMBER OF COLUMNS OF THE MATRIX A. * N MUST BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * ALPHA - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY A MUST * CONTAIN THE MATRIX OF COEFFICIENTS. * UNCHANGED ON EXIT. * * LDA - INTEGER. * ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDA MUST BE AT LEAST * MAX( 1, M ). * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST * ( 1 + ( N - 1 )*ABS( INCX ) ) WHEN TRANS = 'N' OR 'N' * AND AT LEAST * ( 1 + ( M - 1 )*ABS( INCX ) ) OTHERWISE. * BEFORE ENTRY, THE INCREMENTED ARRAY X MUST CONTAIN THE * VECTOR X. * UNCHANGED ON EXIT. * * INCX - INTEGER. * ON ENTRY, INCX SPECIFIES THE INCREMENT FOR THE ELEMENTS OF * X. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * BETA - DOUBLE PRECISION. * ON ENTRY, BETA SPECIFIES THE SCALAR BETA. WHEN BETA IS * SUPPLIED AS ZERO THEN Y NEED NOT BE SET ON INPUT. * UNCHANGED ON EXIT. * * Y - DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST * ( 1 + ( M - 1 )*ABS( INCY ) ) WHEN TRANS = 'N' OR 'N' * AND AT LEAST * ( 1 + ( N - 1 )*ABS( INCY ) ) OTHERWISE. * BEFORE ENTRY WITH BETA NON-ZERO, THE INCREMENTED ARRAY Y * MUST CONTAIN THE VECTOR Y. ON EXIT, Y IS OVERWRITTEN BY THE * UPDATED VECTOR Y. * * INCY - INTEGER. * ON ENTRY, INCY SPECIFIES THE INCREMENT FOR THE ELEMENTS OF * Y. INCY MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME 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 C 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 SUBROUTINE TRED2(NM,N,A,D,E,Z) C INTEGER I,J,K,L,N,II,NM,JP1 DOUBLE PRECISION A(NM,N),D(N),E(N),Z(NM,N) DOUBLE PRECISION F,G,H,HH,SCALE C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A C SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING C ORTHOGONAL SIMILARITY TRANSFORMATIONS. 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 A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE C LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. C C ON OUTPUT C C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. C C Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX C PRODUCED IN THE REDUCTION. C C A AND Z MAY COINCIDE. IF DISTINCT, A IS UNALTERED. 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 DO 100 I = 1, N C DO 80 J = I, N 80 Z(J,I) = A(J,I) C D(I) = A(N,I) 100 CONTINUE C IF (N .EQ. 1) GO TO 510 C .......... FOR I=N STEP -1 UNTIL 2 DO -- .......... DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... DO 120 K = 1, L 120 SCALE = SCALE + DABS(D(K)) C IF (SCALE .NE. 0.0D0) GO TO 140 130 E(I) = D(L) C DO 135 J = 1, L D(J) = Z(L,J) Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 135 CONTINUE C GO TO 290 C 140 DO 150 K = 1, L D(K) = D(K) / SCALE H = H + D(K) * D(K) 150 CONTINUE C F = D(L) G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G D(L) = F - G C .......... FORM A*U .......... DO 170 J = 1, L 170 E(J) = 0.0D0 C DO 240 J = 1, L F = D(J) Z(J,I) = F G = E(J) + Z(J,J) * F JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L G = G + Z(K,J) * D(K) E(K) = E(K) + Z(K,J) * F 200 CONTINUE C 220 E(J) = G 240 CONTINUE C .......... FORM P .......... F = 0.0D0 C DO 245 J = 1, L E(J) = E(J) / H F = F + E(J) * D(J) 245 CONTINUE C HH = F / (H + H) C .......... FORM Q .......... DO 250 J = 1, L 250 E(J) = E(J) - HH * D(J) C .......... FORM REDUCED A .......... DO 280 J = 1, L F = D(J) G = E(J) C DO 260 K = J, L 260 Z(K,J) = Z(K,J) - F * E(K) - G * D(K) C D(J) = Z(L,J) Z(I,J) = 0.0D0 280 CONTINUE C 290 D(I) = H 300 CONTINUE C .......... ACCUMULATION OF TRANSFORMATION MATRICES .......... DO 500 I = 2, N L = I - 1 Z(N,L) = Z(L,L) Z(L,L) = 1.0D0 H = D(I) IF (H .EQ. 0.0D0) GO TO 380 C DO 330 K = 1, L 330 D(K) = Z(K,I) / H C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(K,I) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * D(K) 360 CONTINUE C 380 DO 400 K = 1, L 400 Z(K,I) = 0.0D0 C 500 CONTINUE C 510 DO 520 I = 1, N D(I) = Z(N,I) Z(N,I) = 0.0D0 520 CONTINUE C Z(N,N) = 1.0D0 E(1) = 0.0D0 RETURN END C SUBROUTINE TQL2(NM,N,D,E,Z,IERR) C INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR DOUBLE PRECISION D(N),E(N),Z(NM,N) DOUBLE PRECISION 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 C DOUBLE PRECISION FUNCTION PYTHAG(A,B) DOUBLE PRECISION A,B C C FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW C DOUBLE PRECISION P,R,S,T,U 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 C SUBROUTINE ORTHG(N,F,P,B,LDB,X,LDX,TMP) C---------------------------------------------------------------------- C C.....BLAS2 GRAM-SCHMIDT ORTHOGONALIZATION PROCEDURE 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---------------------------------------------------------------------- INTEGER N,F,P,LDB,LDX DOUBLE PRECISION B(LDB,1), X(LDX,1),TMP(1), DDOT, S, T INTEGER FP1,FPP, K, KM1 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 C DOUBLE PRECISION 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 DOUBLE PRECISION 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 C SUBROUTINE OPAT (X,Y) C C.....MULTIPLICATION OF TRANSPOSE (NROW BY NCOL SPARSE MATRIX A) BY THE C VECTOR X, STORE IN Y PARAMETER(NMAX=500,NZMAX=2000) DOUBLE PRECISION VALUE(NZMAX), ALPHA DOUBLE PRECISION X(1), Y(1) INTEGER POINTR(NMAX), ROWIND(NZMAX) INTEGER I, J, NCOL, NROW, NMAX, NZMAX COMMON /SPARSE/ VALUE,ALPHA,POINTR,ROWIND,NCOL,NROW C C--------------------- C DO 55 I=1,NCOL Y(I)=0.0D0 55 CONTINUE C DO 10 I =1,NCOL C DO 20 J=POINTR(I),POINTR(I+1)-1 Y(I)=Y(I)+VALUE(J)*X(ROWIND(J)) 20 CONTINUE C 10 CONTINUE RETURN END C SUBROUTINE OPB(N, X, Y, SHIFT) C C.....MULTIPLICATION OF A SHIFTED 2-CYCLIC MATRIX C BY A VECTOR X , WHERE C C C = [ D A ] C [ A' D ] , WHERE A IS NROW BY NCOL (NROW >> NCOL), C AND D = (ALPHA-SHIFT)*I, ALPHA AN UPPER BOUND C FOR THE LARGEST SINGULAR VALUE OF A, AND C SHIFT IS AN APPROXIMATE SINGULAR VALUE OF A. C C HENCE, C IS OF ORDER N:= NROW+NCOL (Y STORES PRODUCT VECTOR) C PARAMETER(NMAX=500,NZMAX=2000) DOUBLE PRECISION VALUE(NZMAX), ALPHA DOUBLE PRECISION X(1), Y(1), SHIFT INTEGER POINTR(NMAX), ROWIND(NZMAX) INTEGER I, J, NROW, NCOL, NMAX, NZMAX COMMON /SPARSE/ VALUE,ALPHA,POINTR,ROWIND,NCOL,NROW INTEGER N C C--------------------- C DO 55 I=1,N Y(I)=(ALPHA-SHIFT)*X(I) 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(NROW+I) 10 CONTINUE C 15 CONTINUE C DO 25 I =NROW+1,N C DO 20 J=POINTR(I-NROW),POINTR(I-NROW+1)-1 Y(I)=Y(I)+VALUE(J)*X(ROWIND(J)) 20 CONTINUE 25 CONTINUE C RETURN END SUBROUTINE DISK(N,SIG,RAD,CSIZE,CLUS) C C.... MONITOR SEPARATION OF GERSHGORIN DISKS C DOUBLE PRECISION SIG(N),RAD(N),RADI,RADIP1,TMP INTEGER CPTR,CSIZE(N),CLUS(N),N,I,K C.... C.... ASSUME ORDERING: SIG(1)<=SIG(2)<=.....<=SIG(N) FOR SIG ARRAY IN TSVD1 C.... C.... RAD:= RADII FOR APPROXIMATE E-VALUES C.... DO 10 I = 1,N-1 RADI=RAD(I) RADIP1=RAD(I+1) TMP=(RADI+RADIP1) - (SIG(I+1)-SIG(I)) C IF(TMP.LE.0.0D0) THEN CLUS(I)=1 ELSE CLUS(I)=0 ENDIF C 10 CONTINUE C C.....CLUSTERING VECTOR FILLED, NOW LOCATE CLUSTERS AND THEIR SIZE C DO 5 I = 1,N CSIZE(I)=1 5 CONTINUE C DO 20 I = 1,N-1 CPTR=I DO 25 K = I-1,1,-1 IF(CSIZE(K).NE.0) THEN GO TO 26 ELSE CPTR=K ENDIF 25 CONTINUE C 26 IF(CLUS(I).EQ.0) THEN IF(CSIZE(I).NE.0) THEN CSIZE(I) = CSIZE(I) +1 ELSE CSIZE(CPTR-1) = CSIZE(CPTR-1) +1 ENDIF C CSIZE(I+1)=CSIZE(I+1)-1 ENDIF C 20 CONTINUE C C.....CSIZE ARRAY INDICATES LOCATION AND SIZE OF CLUSTERS: C C CSIZE(I)=K (K.NE.0) := CLUSTER OF SIZE K WITH FIRST DISK C CENTERED AT SIG(I). C RETURN END SUBROUTINE ISOL(I,RESID,TOL,INIT, * IRES0,CRES0,S) INTEGER S,I DOUBLE PRECISION IRES0(S), CRES0(S), RESID, TOL LOGICAL INIT(S) C.... C.... MONITOR RESIDUAL REDUCTION FOR ISOLATED SINGULAR VALUE APPROXIMATIONS C.... C.... ASSUME ORDERING: SIG(1)<=SIG(2)<=.....<=SIG(S) FOR SIG ARRAY IN TSVD1 C.... C.... RESID:= 2-NORM OF B-EIGENPAIR RESIDUAL C.... IF(IRES0(I).LT.0.0D0.OR.CRES0(I).GE.0.0D0) THEN C IF(CRES0(I).LT.0.0D0) THEN IRES0(I)=RESID ELSE IRES0(I)=RESID CRES0(I)=-1.0D0 ENDIF C ELSE C IF(RESID.LE.TOL*IRES0(I)) INIT(I)=.TRUE. C ENDIF RETURN END SUBROUTINE CLUS(I,SIZE,RESID,TOL,INIT,IRES0,CRES0,TMP,S) INTEGER SIZE,S,I,K DOUBLE PRECISION RESID(SIZE),TMP(SIZE),CRES0(S),IRES0(S) DOUBLE PRECISION DASUM,ERROR,TOL LOGICAL INIT(S) C.... C.... MONITOR RESIDUAL REDUCTION FOR CLUSTERED SINGULAR VALUE APPROXIMATIONS C.... C.... ASSUME ORDERING: SIG(1)<=SIG(2)<=.....<=SIG(N) FOR SIG ARRAY IN TSVD1 C.... C.... RESID:= 2-NORM OF B-EIGENPAIR RESIDUALS C.... I:= FIRST DISK IN CLUSTER C.... SIZE:= NUMBER OF DISKS IN CLUSTER C.... DO 20 K=1,SIZE TMP(K)=RESID(K)*RESID(K) 20 CONTINUE C ERROR=DSQRT(DASUM(SIZE,TMP,1)) C IF(CRES0(I).LT.0.0D0) THEN IF(IRES0(I).LT.0.0D0) THEN CRES0(I)=ERROR ELSE CRES0(I)=ERROR IRES0(I)=-1.0D0 ENDIF C ELSE C IF(ERROR.LE.TOL*CRES0(I)) THEN DO 40 K=I,I+SIZE-1 INIT(K)=.TRUE. 40 CONTINUE ENDIF ENDIF RETURN END SUBROUTINE CGT(N,LEFT,W,V,LDV,CGITER,SIG,SIGMAX, * KOUNT,EPS) C C.... CG FOR INDEPENDENT SYSTEMS IN TRACE MIN. OPTIMIZATION STEP C.... C.... V STORES CURRENT ORTHOG BASIS FOR R(Y) C.... SET NCG TO BE AT LEAST N = ORDER OF MATRIX B C.... PARAMETER(NCG=500) DOUBLE PRECISION W(N), V(LDV,LEFT), * Z(NCG),V2(NCG),R(NCG),R2(NCG), * P(NCG),Z2(NCG),SIG,SIGMAX,EPS, * DENOM,DABS,DDOT,BOUND,A,A0, * B,DNRM2,ERROR,RNORM,RNORM0,VNORM INTEGER CGITER,N,MAXI,NCG,LEFT,LDV,KOUNT,I,II,J,K C MAXI=N KOUNT=0 BOUND=(SIG/SIGMAX)**2 C C.... W0=W BY DEFINITION, GET FIRST RESIDUAL VIA ORTHOG. PROJECTOR C CALL OPB(N,W,Z,0.0D0) KOUNT=KOUNT+1 CALL DGEMV('T',N,LEFT,1.0D0,V,LDV,Z,1,0.0D0,R,1) CALL DGEMV('N',N,LEFT,1.0D0,V,LDV,R,1,0.0D0,V2,1) C DO 10 I=1,N R(I)=Z(I)-V2(I) 10 CONTINUE C RNORM0=DNRM2(N,R,1) C DO 20 I=1,N P(I)=R(I) 20 CONTINUE C C.....(MAIN ITERATION LOOP 30) C DO 30 I = 1,MAXI C CALL OPB(N,P,V2,0.0D0) KOUNT=KOUNT+1 C DENOM=DDOT(N,P,1,V2,1) IF(DENOM.LE.0.0) RETURN A = DDOT(N,R,1,R,1)/DENOM IF(I.EQ.1) A0=A CALL DAXPY(N,-A,P,1,W,1) C DO 40 J=1,N R2(J)=R(J) 40 CONTINUE CALL DGEMV('T',N,LEFT,1.0D0,V,LDV,V2,1,0.0D0,Z,1) CALL DGEMV('N',N,LEFT,1.0D0,V,LDV,Z,1,0.0D0,Z2,1) DO 45 II=1,N V2(II)=V2(II)-Z2(II) 45 CONTINUE CALL DAXPY(N,-A,V2,1,R2,1) RNORM=DNRM2(N,R2,1) DO 145 K=1,N V2(K)=P(K) 145 CONTINUE CALL DSCAL(N,A,V2,1) VNORM=DNRM2(N,V2,1) C C.... EARLY TERMINATION CODE: C ERROR=DABS(A*RNORM*RNORM/(A0*RNORM0*RNORM0)) IF(ERROR.LE.BOUND.OR.VNORM.LE.EPS) THEN CGITER=CGITER+I GO TO 999 ELSE IF(I.EQ.MAXI) THEN CGITER=CGITER+MAXI GO TO 998 ENDIF C DO 50 J=1,N V2(J)=R2(J) 50 CONTINUE B = DDOT(N,R2,1,R2,1)/DDOT(N,R,1,R,1) CALL DAXPY(N,B,P,1,V2,1) DO 60 J=1,N P(J)=V2(J) 60 CONTINUE DO 65 J=1,N R(J)=R2(J) 65 CONTINUE 30 CONTINUE 998 CONTINUE WRITE(*,*) 'CGTS FAILED TO CONVERGE IN ',MAXI,' ITERATIONS.' 999 CONTINUE RETURN END SUBROUTINE CGTS(N,LEFT,W,V,LDV,CGITER,SIG,SIGOLD,SIGMAX, * SHIFT,KOUNT,EPS) C C.... CG FOR INDEPENDENT SYSTEMS IN TRACE MIN. OPTIMIZATION STEP C.... (SHIFT INCORPORATED) C.... V STORES CURRENT ORTHOG BASIS FOR R(Y) C.... C.... SET NCG TO BE AT LEAST N = ORDER OF MATRIX B C.... PARAMETER(NCG=500) DOUBLE PRECISION W(N), V(LDV,LEFT), * Z(NCG),V2(NCG),R(NCG),R2(NCG), * P(NCG),Z2(NCG),SIG,SIGMAX,EPS,SHIFT, * DENOM,DABS,DDOT,BOUND,A,A0,SIGOLD, * B,DNRM2,ERROR,RNORM,RNORM0,VNORM INTEGER CGITER,N,MAXI,NCG,LEFT,LDV,KOUNT,I,II,J,K C MAXI=N KOUNT=0 IF(SIG.NE.SHIFT) THEN BOUND=((SIG-SHIFT)/(SIGMAX-SHIFT))**2 ELSE BOUND=((SIGOLD-SHIFT)/(SIGMAX-SHIFT))**2 ENDIF C C.... W0=W BY DEFINITION, GET FIRST RESIDUAL VIA ORTHOG. PROJECTOR C CALL OPB(N,W,Z,SHIFT) KOUNT=KOUNT+1 CALL DGEMV('T',N,LEFT,1.0D0,V,LDV,Z,1,0.0D0,R,1) CALL DGEMV('N',N,LEFT,1.0D0,V,LDV,R,1,0.0D0,V2,1) C DO 10 I=1,N R(I)=Z(I)-V2(I) 10 CONTINUE C RNORM0=DNRM2(N,R,1) C DO 20 I=1,N P(I)=R(I) 20 CONTINUE C C.... (MAIN ITERATION LOOP 30) C DO 30 I = 1,MAXI CALL OPB(N,P,V2,SHIFT) KOUNT=KOUNT+1 DENOM=DDOT(N,P,1,V2,1) IF(DENOM.LE.0.0) RETURN A = DDOT(N,R,1,R,1)/DENOM IF(I.EQ.1) A0=A CALL DAXPY(N,-A,P,1,W,1) C DO 40 J=1,N R2(J)=R(J) 40 CONTINUE CALL DGEMV('T',N,LEFT,1.0D0,V,LDV,V2,1,0.0D0,Z,1) CALL DGEMV('N',N,LEFT,1.0D0,V,LDV,Z,1,0.0D0,Z2,1) DO 45 II=1,N V2(II)=V2(II)-Z2(II) 45 CONTINUE CALL DAXPY(N,-A,V2,1,R2,1) RNORM=DNRM2(N,R2,1) DO 145 K=1,N V2(K)=P(K) 145 CONTINUE CALL DSCAL(N,A,V2,1) VNORM=DNRM2(N,V2,1) C C.... EARLY TERMINATION CODE: C ERROR=DABS(A*RNORM*RNORM/(A0*RNORM0*RNORM0)) IF(ERROR.LE.BOUND.OR.VNORM.LE.EPS) THEN CGITER=CGITER+I GO TO 999 ELSE IF(I.EQ.MAXI) THEN CGITER=CGITER+MAXI GO TO 998 ENDIF C DO 50 J=1,N V2(J)=R2(J) 50 CONTINUE B = DDOT(N,R2,1,R2,1)/DDOT(N,R,1,R,1) CALL DAXPY(N,B,P,1,V2,1) DO 60 J=1,N P(J)=V2(J) 60 CONTINUE DO 65 J=1,N R(J)=R2(J) 65 CONTINUE C 30 CONTINUE 998 CONTINUE WRITE(*,*) 'CGTS FAILED TO CONVERGE IN ',MAXI,' ITERATIONS.' 999 CONTINUE RETURN END .