C-----LSMULT------------------------------------------------------------ C C CONTAINS SUBROUTINES LANCZS, USPECS, STRAN, AND SVMAT C FOR USE WITH THE LANCZOS SINGULAR VALUE/VECTOR PROGRAMS C C NONPORTABLE CONSTRUCTIONS: C 1. THE ENTRY MECHANISM USED TO PASS THE STORAGE LOCATIONS C OF THE USER-SPECIFIED MATRIX FROM THE SUBROUTINE USPEC C TO THE MATRIX-VECTOR MULTIPLY SUBROUTINES SVMAT AND C STRAN. C 2. IN THE SAMPLE USPEC PROVIDED: THE FREE FORMAT (8,*), C AND THE FORMAT (20A4). C C-----START OF LANCZS--------------------------------------------------- C SUBROUTINE LANCZS(MATVEC,MTRAN,BETA,V1,V2,G,KMAX,MOLD1, 1 M,N,IPAR,IIX) C C----------------------------------------------------------------------- DOUBLE PRECISION BETA(1),V1(1),V2(1),SUM,TEMP,ONE,ZERO REAL G(1) DOUBLE PRECISION FINPRO INTEGER IPAR EXTERNAL MATVEC,MTRAN C----------------------------------------------------------------------- C COMPUTE T(1,MEV) FOR SYMMETRIZED VERSION OF GIVEN A-MATRIX. C C ---- ---- C | 0 A | C B = | | C | A-TRANSPOSE 0 | C ---- ---- C C WHERE A IS AN M BY N REAL SPARSE MATRIX, USING STARTING C VECTORS OF THE FORM (V1,0) WHEN THE FLAG IPAR = 2 AND C OF THE FORM (0,V2) WHEN THE FLAG IPAR = 1. V1 IS OF C DIMENSION M, THE ROW DIMENSION OF A, AND V2 IS OF DIMENSION C N, THE COLUMN DIMENSION OF A. C C WITH STARTING VECTORS OF THESE FORMS, THE LANCZOS VECTORS C GENERATED ALTERNATE BETWEEN THESE 2 FORMS AND ALL OF THE C DIAGONAL ENTRIES OF THE LANCZOS TRIDIAGONAL MATRICES T(1,MEV) C GENERATED ARE 0. C C LANCZS USES 2 USER-SUPPLIED SUBROUTINES MATVEC AND MTRAN. C MAIN PROGRAM CALLS THESE SVMAT AND STRAN, RESPECTIVELY. C CALLING SEQUENCES ARE C C CALL MATVEC(V2,V1,SUM) C CALL MTRAN(V1,V2,SUM) C C MATVEC COMPUTES V1 = A*V2 - SUM*V1. C MTRAN COMPUTES V2 = (A-TRANSPOSE)*V1 - SUM*V2. C C ON EXIT V1 AND V2 CONTAIN THE NONZERO PARTS OF THE C LAST TWO LANCZOS VECTORS. C C IF MOLD1 = 1 THEN T(1,KMAX) IS GENERATED FROM SCRATCH. C IF MOLD1 > 1 THEN A PREVIOUSLY-GENERATED T-MATRIX OF SIZE C (MOLD1-1) IS EXTENDED TO ONE OF SIZE KMAX. SINGULAR VALUE C PRGORAMS CAN ONLY UTILIZE T-MATRICES OF EVEN ORDER. C BETA(KMAX+1) IS ALSO COMPUTED FOR USE IN THE ERROR ESTIMATES. C C----------------------------------------------------------------------- ONE = 1.0D0 ZERO = 0.0D0 ITNUM = MOLD1 C IF (ITNUM .GT. 1) GO TO (80,100), IPAR C C NO PREVIOUS BETA HISTORY BETA(1) = ZERO IIL = IIX IF (IPAR .EQ. 2) GO TO 40 C C----------------------------------------------------------------------- C IPAR = 1 SO SET V2 EQUAL TO A UNIT RANDOM VECTOR AND SET V1 = 0. CALL GENRAN(IIL,G,N) C----------------------------------------------------------------------- C DO 10 J = 1,N 10 V2(J) = G(J) C C----------------------------------------------------------------------- TEMP = FINPRO(N,V2(1),1,V2(1),1) C----------------------------------------------------------------------- C SUM = ONE/DSQRT(TEMP) DO 20 J = 1,M 20 V1(J) = ZERO C DO 30 J = 1,N 30 V2(J) = V2(J)*SUM GO TO 100 C 40 CONTINUE C C----------------------------------------------------------------------- C IPAR = 2 SO SET V1 EQUAL TO A UNIT RANDOM VECTOR AND SET V2 = 0. CALL GENRAN(IIL,G,M) C----------------------------------------------------------------------- C DO 50 J=1,M 50 V1(J) = G(J) C C----------------------------------------------------------------------- TEMP = FINPRO(M,V1(1),1,V1(1),1) C----------------------------------------------------------------------- C SUM = ONE/DSQRT(TEMP) DO 60 J = 1,N 60 V2(J) = ZERO DO 70 J = 1,M 70 V1(J) = V1(J)*SUM C C BELOW IS START FOR MOLD1 > 1 AND IPAR = 1 C DO ONE ITERATION OF LANCZOS TO OBTAIN (0,V2) C 80 CONTINUE SUM = BETA(ITNUM) C C----------------------------------------------------------------------- CALL MTRAN(V1,V2,SUM) C----------------------------------------------------------------------- C C----------------------------------------------------------------------- SUM = FINPRO(N,V2(1),1,V2(1),1) C----------------------------------------------------------------------- C ITNUM = ITNUM + 1 BETA(ITNUM) = DSQRT(SUM) SUM = ONE/BETA(ITNUM) C DO 90 J = 1,N 90 V2(J) = V2(J)*SUM C IPAR = 2 IF (ITNUM .GT. KMAX) GO TO 120 C C BELOW IS START FOR MOLD1 > 1 AND IPAR = 2 C DO ONE ITERATION OF LANCZOS TO OBTAIN (V1,0) C 100 CONTINUE SUM = BETA(ITNUM) C C----------------------------------------------------------------------- CALL MATVEC(V2,V1,SUM) C----------------------------------------------------------------------- C C----------------------------------------------------------------------- SUM = FINPRO(M,V1(1),1,V1(1),1) C----------------------------------------------------------------------- C ITNUM = ITNUM + 1 BETA(ITNUM) = DSQRT(SUM) SUM = ONE/BETA(ITNUM) C DO 110 J = 1,M 110 V1(J)= V1(J) * SUM C IPAR = 1 IF (ITNUM .GT. KMAX) GO TO 120 GO TO 80 C 120 CONTINUE C RETURN C-----END OF LANCZS----------------------------------------------------- END C C-----START OF USPEC (GENERAL SPARSE, RECTANGULAR MATRIX)--------------- C C SUBROUTINE USPEC(M,N,MATNO) SUBROUTINE SUSPEC(M,N,MATNO) C C-------------------------------------------------------------------- DOUBLE PRECISION A(10000) INTEGER IROW(10000),ICOL(3010) C-------------------------------------------------------------------- C DIMENSIONS ARRAYS NEEDED TO DEFINE THE USER-SUPPLIED C M X N RECTANGULAR A-MATRIX, READS IN VALUES OF THESE C ARRAYS AND THEN PASSES THE STORAGE LOCATIONS OF THESE C ARRAYS TO THE CORRESPONDING MATRIX-VECTOR MULTIPLY C SUBROUTINES SVMAT AND STRAN. C C THE A-MATRIX IS STORED IN THE FOLLOWING SPARSE FORMAT: C M = NUMBER OF ROWS IN A. C N = NUMBER OF COLUMNS IN A. C NZ = NUMBER OF NONZERO ENTRIES IN A-MATRIX. C ICOL(J), J=1,N IS NUMBER OF NONZERO ENTRIES IN COLUMN J. C IROW(K), K = 1,NZ IS THE ROW INDEX FOR CORRESPONDING A(K). C A(K), K=1,NZ IS NONZERO ENTRIES IN A, COLUMN BY COLUMN. C IT IS ASSUMED THAT ICOL(J) > 0 FOR ALL J C C NOTE: ASSOCIATED SUBROUTINES SVMAT AND STRAN ASSUME THAT C M >= N. C C----------------------------------------------------------------------- C READ IN MATRIX FROM FILE 8 C READ(8,10) NZ,MOLD,NOLD,MATOLD 10 FORMAT(I10,2I6,I8) C WRITE(6,20) NZ,MOLD,NOLD,MATOLD 20 FORMAT(6X,'NZ',4X,'MOLD',4X,'NOLD',4X,'MATOLD'/I10,2I6,I10/) C C TEST OF PARAMETER CORRECTNESS ITEMP = (MOLD-M)**2 + (NOLD-N)**2 + (MATOLD-MATNO)**2 C IF (ITEMP.EQ.0) GO TO 40 C WRITE(6,30) 30 FORMAT(' PROGRAM TERMINATES BECAUSE EITHER ORDERS OF OR LABELS FOR 1 MATRIX DISAGREE') GO TO 70 C 40 CONTINUE C C NUMBER OF NONZERO ENTRIES IN EACH COLUMN IS READ IN C THEN THE CORRESPONDING ROW INDEX FOR EACH SUCH ENTRY IS READ READ(8,50) (ICOL(K), K=1,N) READ(8,50) (IROW(K), K=1,NZ) 50 FORMAT(13I6) C C READ IN THE NONZERO ENTRIES IN THE MATRIX READ(8,60) (A(K), K=1,NZ) 60 FORMAT(3E25.16) C 50 FORMAT(4E19.10) C C----------------------------------------------------------------------- C PASS STORAGE LOCATIONS OF ARRAYS THAT DEFINE THE MATRIX TO C THE MATRIX-VECTOR MULTIPLY SUBROUTINES SVMAT AND STRAN CALL SMATVE(A,ICOL,IROW,M,N) CALL STRANE(A,ICOL,IROW,M,N) C----------------------------------------------------------------------- C C-----END OF USPEC------------------------------------------------------ RETURN 70 STOP END C C-----STRAN (GENERAL SPARSE MATRIX)------------------------------------- C C SUBROUTINE STRAN(W,U,SUM) SUBROUTINE SSTRAN(W,U,SUM) C C----------------------------------------------------------------------- DOUBLE PRECISION W(1),U(1),A(1),SUM,TEMP INTEGER IROW(1),ICOL(1) C----------------------------------------------------------------------- C SUBROUTINE TO COMPUTE U = (A-TRANSPOSE)*W - SUM*U WHERE A IS C A GENERAL, SPARSE M X N MATRIX WITH M >= N. C C ASSUMES MATRIX IS STORED IN SPARSE FORMAT GIVEN IN C CORRESPONDING USPEC SUBROUTINE. C----------------------------------------------------------------------- JLAST = 0 DO 20 J = 1,N JFIRST = JLAST + 1 JLAST = JLAST + ICOL(J) TEMP = -SUM*U(J) C DO 10 K = JFIRST,JLAST IK = IROW(K) 10 TEMP = A(K)*W(IK) + TEMP C 20 U(J) = TEMP C RETURN C C----------------------------------------------------------------------- ENTRY STRANE(A,ICOL,IROW,M,N) C----------------------------------------------------------------------- C C-----END OF STRAN FOR GENERAL SPARSE MATRIX---------------------------- RETURN END C C-----SVMAT (GENERAL SPARSE MATRIX)------------------------------------- C C SUBROUTINE SVMAT(W,U,SUM) SUBROUTINE SSVMAT(W,U,SUM) C C----------------------------------------------------------------------- DOUBLE PRECISION W(1),U(1),A(1),SUM,TEMP INTEGER IROW(1),ICOL(1) C----------------------------------------------------------------------- C SUBROUTINE TO COMPUTE U = A*W - SUM*U WHERE A IS A C GENERAL, SPARSE M X N MATRIX WITH M >= N. C C ASSUMES THAT THE MATRIX IS STORED IN THE SPARSE FORMAT C GIVEN IN THE CORRESPONDING USPEC SUBROUTINE. C----------------------------------------------------------------------- DO 10 I = 1,M 10 U(I) = -SUM*U(I) C C MAIN LOOP. PROCESSING PROCEEDS COL BY COL. JFIRST AND JLAST ARE C POINTERS TO THE FIRST AND LAST NONZEROS IN COLUMN J. C JLAST = 0 DO 30 J = 1,N JFIRST = JLAST + 1 JLAST = JLAST + ICOL(J) TEMP = W(J) C DO 20 K = JFIRST,JLAST IK = IROW(K) 20 U(IK) = U(IK) + A(K)*TEMP C 30 CONTINUE C RETURN C C----------------------------------------------------------------------- ENTRY SMATVE(A,ICOL,IROW,M,N) C----------------------------------------------------------------------- C C----END OF SVMAT FOR GENERAL SPARSE MATRICES--------------------------- RETURN END C C-----ROUTINES FOR 'DIAGONAL' TEST MATRICES----------------------------- C DMATV,DMTRAN,DIAGSP SUBROUTINES ARE FOR RECTANGULAR DIAGONAL C TEST MATRICES. C C-----START OF USPEC FOR 'DIAGONAL' TEST MATRIX------------------------- C SUBROUTINE USPEC(M,N,MATNO) C SUBROUTINE DIAGSP(M,N,MATNO) C C DEFINES 'DIAGONAL' MATRIX OF FOLLOWING FORM C C ----- ----- C | 0 0 D | C A = | 0 0 0 | C |D-TRANS 0 0 | C ----- ----- C C WHERE D IS DIAGONAL MATRIX OF ORDER N, AND IN THE C MIDDLE THERE ARE (M-N) ROWS OF ZEROES. C CALLS ENTRY TO MATRIX-VECTOR MULTIPLY SUBROUTINE TO PASS C STORAGE LOCATION OF THE D-ARRAY AND THE ORDERS M AND N. C C NOTE: ASSOCIATED MATRIX-VECTOR SUBROUTINES ASSUME THAT C M >= N. C----------------------------------------------------------------------- DOUBLE PRECISION D(1000), SPACE REAL EXPLAN(20) C----------------------------------------------------------------------- C READ(8,10) EXPLAN 10 FORMAT(20A4) READ(8,*) MOLD,NOLD,NUNIF,SPACE,D(1) C IF(N.NE.NOLD.OR.M.NE.MOLD) GO TO 80 C COMPUTE THE UNIFORM PORTION OF THE SPECTRUM DO 20 J=2,NUNIF 20 D(J) = D(1) - DFLOAT(J-1)*SPACE NUNIF1=NUNIF + 1 READ(8,10) EXPLAN DO 30 J=NUNIF1,N 30 READ(8,*) D(J) NNUNIF = NOLD - NUNIF WRITE(6,40) NOLD,SPACE,NNUNIF,D(1) 40 FORMAT(/' DIAGONAL TEST MATRIX, SIZE = ',I4/' MOST ENTRIES ARE ', 1E10.3,' UNITS APART.',I3,' ENTRIES'/' ARE IRREGULARLY SPACED. FIRS 1T ENTRY IS ',E10.3/) NB = NUNIF - 2 C C C PRINT OUT DIAGONAL PORTION OF A-MATRIX WRITE(6,50) (D(I), I=1,10 ) WRITE(6,60) (D(I), I = NB,N) MNDIF = MOLD - NOLD IF(MNDIF.NE.0) WRITE(6,70) MNDIF 50 FORMAT(/' SINGULAR VALUE LANCZOS TEST, 1ST 10 ENTRIES OF DIAGONAL 1A-MATRIX = '/(3E22.14)) 60 FORMAT(/' MIDDLE UNIFORM PORTION OF MATRIX IS NOT PRINTED OUT'/ 1' END OF UNIFORM PLUS NONUNIFORM SECTION = '/(3E22.14)) 70 FORMAT(I4,' ZERO ROWS ARE ADDED TO THE DIAGONAL TO MAKE IT RECTANG 1ULAR'/) C C DIAGONAL GENERATION COMPLETE C C----------------------------------------------------------------------- C CALL ENTRY TO MATRIX-VECTOR MULTIPLY SUBROUTINES TO PASS C STORAGE LOCATION OF D-ARRAY AND ORDER OF A-MATRIX. CALL DMATVE(D,M,N) CALL DMTRAE(D,M,N) C----------------------------------------------------------------------- C RETURN 80 WRITE(6,90) MOLD,NOLD,M,N 90 FORMAT(' PROGRAM TERMINATES MOLD=',I5,' N.E. M=',I5,' OR NOLD=', 1I5,' N.E. N=',I5) C-----END OF USPEC SUBROUTINE FOR 'DIAGONAL' TEST MATRICES-------------- STOP END C C-----DSVMAT ('DIAGONAL' TEST MATRICES)--------------------------------- C C SUBROUTINE DSVMAT(Z,W,SUM) SUBROUTINE SVMAT(Z,W,SUM) C C----------------------------------------------------------------------- DOUBLE PRECISION A(1),Z(1),W(1),SUM C----------------------------------------------------------------------- C C COMPUTES W = A*Z - SUM*W . ASSUMES THAT M >= N. DO 10 I = 1,N 10 W(I) = A(I)*Z(I) - SUM *W(I) IF(M.EQ.N) RETURN N1 = N+1 DO 20 I = N1,M 20 W(I) = -SUM*W(I) RETURN C C----------------------------------------------------------------------- C STORAGE LOCATIONS OF THE A-ARRAY C AND THE ORDER OF THE A-MATRIX ARE PASSED TO THE MATVEC SUBROUTINE. C ENTRY MATVE(A,M,N) ENTRY DMATVE(A,M,N) C----------------------------------------------------------------------- C C-----END OF MATRIX -VECTOR MULTIPLY 'DIAGONAL' TEST PROBLEMS----------- RETURN END C C-----MATRIX-VECTOR MULTIPLY FOR 'DIAGONAL' TEST MATRICES--------------- C SUBROUTINE STRAN(Z,W,SUM) C SUBROUTINE DSTRAN(Z,W,SUM) C C----------------------------------------------------------------------- DOUBLE PRECISION A(1),Z(1),W(1),SUM C----------------------------------------------------------------------- C C COMPUTES W = A-TRANSPOSE*Z - SUM*W . ASSUMES M >= N. DO 10 I = 1,N 10 W(I) = A(I)*Z(I)- SUM*W(I) RETURN C C----------------------------------------------------------------------- C STORAGE LOCATIONS OF THE A-ARRAY AND THE ORDER C OF THE A-MATRIX ARE OBTAINED FROM USPEC SUBROUTINE. C ENTRY MTRANE(A,M,N) ENTRY DMTRAE(A,M,N) C----------------------------------------------------------------------- C C-----END OF SPARSE SYMMETRIC MATRIX-VECTOR MULTIPLY-------------------- RETURN END