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 SIS1 C C.... THIS SAMPLE PROGRAM USES RITZIT TO COMPUTE SINGULAR TRIPLETS OF A VIA C.... THE EQUIVALENT SYMMETRIC EIGENVALUE PROBLEM C.... C.... B X = LAMBDA X, WHERE X' = (U',V'), LAMBDA = ALPHA +/- SIGMA, C.... C.... [ALPHA*I A] C.... B = [ ] , AND 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. C.... IN THIS 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 VALUE OF THE MATRIX A. C.... (A' = TRANSPOSE OF A) C.... C.... USER SUPPLIED ROUTINES: OPB, OPC C.... C.... OPB(NROW+NCOL,X,Y) TAKES AN (M+N)-VECTOR X AND SHOULD RETURN B*X IN Y. C.... OPC(NROW+NCOL,X,Y) TAKES AN (M+N)-VECTOR X AND SHOULD RETURN C*X IN Y, C.... WHERE C=[B-ALPHA*I], AND I IS THE (M+N) ORDER IDENTITY MATRIX. C.... C.... USER SHOULD REPLACE CALLS TO ETIME WITH AN APPROPRIATE CALL TO C.... AN INSTRINSIC TIMING ROUTINE THAT RETURNS ELAPSED USER TIME C.... C.... NSIG SPECIFIES MAXIMUM NUMBER OF SINGULAR TRIPLETS DESIRED C.... VECTORS SPECIFIES WHETHER OR NOT VECTORS ARE DESIRED ON OUTPUT C.... PARAMETER(NSIG=100) C.... PARAMETER(NMAX=500,NZMAX=2000) REAL*8 VALUE(NZMAX), ALPHA INTEGER POINTR(NMAX), ROWIND(NZMAX) COMMON /MATRIXA/ NCOL,NROW,VALUE,POINTR,ROWIND,ALPHA C C------------------------------------------------------------ C REAL*8 X(NMAX,NSIG),U(NMAX),W(NMAX,NSIG), 1 D(NSIG),F(NSIG),CX(NMAX),RQ(NSIG),B(NSIG,NSIG), 2 S(NSIG),EPS,XNORM,DABS,DDOT,DSUM,DMAX1, 3 TMP1, TMP2, TMP3, TMP4 C REAL TIME, T0, TMP(2) C C------------------------------------------------------------ C C ----- HARWELL-BOEING COLLECTION DATA STRUCTURE C C------------------------------------------------------------ C INTEGER NROW, NCOL INTEGER NNZERO, NRHS, RESULTS, MATRIX, P, EM, EM2, INF C CHARACTER TITLE*72, KEY*8, TYPE*3, PTRFMT*16, INDFMT*16, 1 VALFMT*20, RHSFMT*20 C CHARACTER*30 NAME LOGICAL VECTORS C C--------------------- C EXTERNAL INTROS,OPB,OPC C COMMON /OUTPUT/ INF COMMON /COUNT/ MXVCOUNT INTEGER MXVCOUNT C INPUT=5 MATRIX=7 RESULTS=9 INF=4 C C OPEN FILES FOR INPUT/OUTPUT C OPEN (MATRIX,FILE='LAI7') OPEN (RESULTS,FILE='SIO9') OPEN (INPUT,FILE='SII5') OPEN (INF,FILE='SIO4') WRITE(INF,9994) 9994 FORMAT('-----------------------------------------------',/ * 'INTERMEDIATE OUTPUT PARMS:',// * 'M := CURRENT DEGREE OF CHEBYSHEV POLYNOMIAL',/ * 'S := NEXT ITERATION STEP ',/ * 'G := NUMBER OF ACCEPTED EIGENVALUES OF B ',/ * 'H := NUMBER OF ACCEPTED EIGENVECTORS OF B ',/ * 'F := VECTOR OF ERRORS FOR EIGENPAIRS OF B ',/ * '-----------------------------------------------',/ * ' M S G H F',/ * '-----------------------------------------------'/) 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(*,*) ' 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 LDA = NMAX N=NROW+NCOL C READ (INPUT,*) NAME,EM,NUMEXTRA,KM,EPS,VECTORS C C NAME := DATASET NAME C EM := NUMBER OF DESIRED TRIPLETS C NUMEXTRA := ADDITIONAL VECTORS TO CARRY C KM := MAXIMUM NUMBER OF ITERATIONS C EPS := TOLERANCE FOR ACCEPTING SINGULAR VECTORS C VECTORS := OUTPUT SINGULAR VECTORS (BOOLEAN) C C P = EM + NUMEXTRA EM2=EM C C.....ESTIMATE ALPHA C DO 75 JJ=1,NCOL DSUM=0.0D0 DO 85 II=POINTR(JJ),POINTR(JJ+1)-1 DSUM=DSUM+VALUE(II) 85 CONTINUE ALPHA=DMAX1(ALPHA,DSUM) 75 CONTINUE C T0 = ETIME(TMP(1)) CALL RITZIT(LDA,N,P,KM,EPS,OPB,INTROS,EM,X,D,U,W,B,F,CX,RQ,S, * IMEM) TIME = ETIME(TMP(1))-T0 C WRITE(RESULTS,2000) N,KM,EM2,EM,P,IMEM,VECTORS,ALPHA,EPS, * TITLE,NAME,NROW,NCOL,N C 9998 FORMAT(/1X,'...... SISVD EXECUTION TIME=',1PE10.2) 9999 FORMAT(1X,'...... ' * /1X,'...... ',16X,' MXV =',I12 * /1X,'...... ' * /1X,'...... ',4X,'COMPUTED SINGULAR VALUES',2X, * '(','RESIDUAL NORMS',')' * /1X,'...... ' * /(1X,'...... ',I3,3X,1PE22.14,2X,'(',1PE11.2,')')) 2000 FORMAT( * 1X,'... ' * /1X,'... SOLVE THE CYCLIC EIGENPROBLEM' * /1X,'... NO. OF EQUATIONS =',I10 * /1X,'... MAX. NO. OF ITERATIONS =',I10 * /1X,'... NO. OF DESIRED EIGENPAIRS =',I10 * /1X,'... NO. OF APPROX. EIGENPAIRS =',I10 * /1X,'... INITIAL SUBSPACE DIM. =',I10 * /1X,'... MEMORY REQUIRED (BYTES) =',I10 * /1X,'... WANT S-VECTORS? [T/F] =',L10 * /1X,'... ALPHA =',1PE10.2 * /1X,'... TOLERANCE =',1PE10.2 * //1X,A50 * / 1X,A50 * /1X,'... NO. OF DOCUMENTS (ROWS) = ',I8 * /1X,'... NO. OF TERMS (COLS) = ',I8 * /1X,'... ORDER OF MATRIX B = ',I8 * /1X,'... ') C IF(VECTORS) THEN IVEC=8 OPEN(IVEC,FILE='SIO8',FORM='UNFORMATTED') ENDIF C T0 = ETIME(TMP(1)) C DO 35 J = 1, EM2 CALL OPC(N,X(1,J),CX) TMP1=DDOT(N,X(1,J),1,CX,1) TMP2=DDOT(NROW,X(1,J),1,X(1,J),1) TMP3=DDOT(NCOL,X(NROW+1,J),1,X(NROW+1,J),1) TMP4=DSQRT(TMP2+TMP3) CALL DAXPY(N,-TMP1,X(1,J),1,CX,1) XNORM=DDOT(N,CX,1,CX,1) XNORM=DSQRT(XNORM) F(J)=XNORM/TMP4 D(J)=DABS(TMP1) IF(VECTORS) WRITE(IVEC) (X(I,J),I=1,N) 35 CONTINUE C TIME = TIME + ETIME(TMP(1))-T0 C WRITE(RESULTS,9998) TIME WRITE(RESULTS,9999)MXVCOUNT,(I,D(I),F(I),I = 1,EM2) C STOP END C SUBROUTINE INTROS(KS,KG,KH,F,M) C C INTROS IS A SUBROUTINE USED TO OBTAIN C INFORMATION OR EXERT CONTROL DURING EXECUTION. INTROS IS CALLED C WITH PARAMETERS (KS,KG,KH,F,M), WHERE KS IS THE NUMBER C OF THE NEXT ITERATION STEP, KG IS THE NUMBER OF ALREADY C ACCEPTED EIGENVECTORS, KH IS THE NUMBER OF ALREADY ACCEPTED C EIGENVALUES, AND F IS THE ARRAY OF ERROR QUANTITIES FOR C VECTORS OF X. AN ELEMENT OF F HAS THE VALUE 4.0 C UNTIL THE CORRESPONDING EIGENVALUE HAS BEEN ACCEPTED. C M IS THE DEGREE OF THE CURRENT CHEBYSHEV POLYNOMIAL. C INTEGER KG,KH,INF,M,KS,L REAL*8 F(1) COMMON /OUTPUT/ INF L = KH + 1 WRITE(INF,10) M,KS,KG,KH,F(1) IF(L.GE.2) WRITE(INF,11) (F(I),I=2,L) 10 FORMAT(/4I8,E15.6) 11 FORMAT( 32X,E15.6) RETURN END C SUBROUTINE RITZIT(NM,N,KP,KM,EPS,OPB,INF,KEM,X,D,U,W,B,F,CX, * RQ,S,IMEM) C INTEGER I,J,K,L,M,N,IG,II,IK,IP,KG,KH,KM,KP,KS,KZ,L1,M1,NM,KEM, * KZ1,KZ2,JP,IMEM REAL*8 X(NM,KP),D(KP),U(NM),W(NM,KP),B(KP,KP),F(KP),CX(KP), * RQ(KP) REAL*8 E,S(KP),T,E1,E2,XK,EE2,EPS,XKM,XKS,XM1 REAL*8 DABS,DEXP,DFLOAT,DLOG,DMAX1,DSQRT,DDOT C INTEGER IABS,IDINT,MAX0,MOD C C CALLS SUBROUTINES OPB,INF,TRED2,TQL2,DGEMM,DGEMV,DAXPY C CALLS FUNCTION DDOT,PYTHAG,LSAME C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE RITZIT, C NUM. MATH. 16, 205-223(1970) BY RUTISHAUSER. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 284-302(1971). C C THIS SUBROUTINE DETERMINES THE ABSOLUTELY LARGEST EIGENVALUES C AND CORRESPONDING EIGENVECTORS OF A REAL SYMMETRIC MATRIX BY C SIMULTANEOUS ITERATION. C C ON INPUT: C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETER X AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT; C C N IS THE ORDER OF THE MATRIX A (MATRIX B FOR SVD PROBLEM); C C KP IS THE NUMBER OF SIMULTANEOUS ITERATION VECTORS; C C KM IS THE MAXIMUM NUMBER OF ITERATION STEPS TO BE PERFORMED. C IF STARTING VALUES FOR THE ITERATION VECTORS ARE AVAILABLE, C KM SHOULD BE PREFIXED WITH A MINUS SIGN; C C EPS IS THE TOLERANCE FOR ACCEPTING EIGENVECTORS; C C OPB IS THE NAME OF THE SUBROUTINE THAT DEFINES THE MATRIX B. C OPB IS CALLED WITH PARAMETERS (N,U,W) AND MUST COMPUTE C W=BU WITHOUT ALTERING U FOR SVD PROBLEM; C C INF IS THE NAME OF THE SUBROUTINE THAT MAY BE USED TO OBTAIN C INFORMATION OR EXERT CONTROL DURING EXECUTION. INF IS CALLED C WITH PARAMETERS (KS,KG,KH,F,KM,KEM), WHERE KS IS THE NUMBER C OF THE NEXT ITERATION STEP, KG IS THE NUMBER OF ALREADY C ACCEPTED EIGENVECTORS, KH IS THE NUMBER OF ALREADY ACCEPTED C EIGENVALUES, AND F IS THE ARRAY OF ERROR QUANTITIES FOR C THE VECTORS OF X. AN ELEMENT OF F HAS THE VALUE 4.0 C UNTIL THE CORRESPONDING EIGENVALUE HAS BEEN ACCEPTED; C C KEM IS THE NUMBER OF EIGENVALUES AND CORRESPONDING C EIGENVECTORS DESIRED. KEM MUST BE LESS THAN KP; C C X CONTAINS, IF KM IS NEGATIVE, THE STARTING VALUES FOR C THE ITERATION VECTORS. C C ON OUTPUT: C C KM IS RESET TO THE MAGNITUDE OF ITS INPUT VALUE; C C KEM IS RESET TO THE NUMBER OF EIGENVALUES AND EIGENVECTORS C ACTUALLY ACCEPTED WITHIN THE LIMIT OF KM STEPS; C C IMEM IS SET TO THE APPROXIMATE NUMBER OF BYTES NEEDED FOR C THIS INVOCATION. C C X CONTAINS IN ITS FIRST KEM COLUMNS ORTHONORMAL EIGENVECTORS C OF A CORRESPONDING TO THE EIGENVALUES IN D. THE REMAINING C COLUMNS CONTAIN APPROXIMATIONS TO FURTHER EIGENVECTORS; C C D CONTAINS IN ITS FIRST KEM POSITIONS THE ABSOLUTELY LARGEST C EIGENVALUES OF A (B). THE REMAINING POSITIONS CONTAIN C APPROXIMATIONS TO SMALLER EIGENVALUES; C C U, W, B, F, CX, S, AND RQ ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY. C MODERNIZATIONS BY M. W. BERRY, UNIV. OF TENNESSEE, 1991. C C ------------------------------------------------------------------ C C.....GET AUXILLARY MEMORY COUNT (8-BYTE WORDS) C (ASSUME NM=N FOR COUNTS) C IMEM = N*KP + KP + KP IMEM = IMEM + N + N*KP + KP*KP + KP + KP + KP IMEM = 8*IMEM C EE2 = 1.0D0 + 1.0D-1 * EPS E = 0.0D0 KG = 0 KH = 0 KZ = 1367 KZ1 = 0 KZ2 = 0 KS = 0 M = 1 C DO 30 L = 1, KP F(L) = 4.0D0 CX(L) = 0.0D0 RQ(L) = 0.0D0 30 CONTINUE C IF (KM .LT. 0) GO TO 60 C :::::::::: GENERATE RANDOM INITIAL ITERATION VECTORS :::::::::: DO 50 L = 1, KP C DO 40 J = 1, N KZ = MOD(3125*KZ,65536) X(J,L) = DFLOAT(KZ-32768) 40 CONTINUE C 50 CONTINUE C 60 KM = IABS(KM) L1 = 1 II = 1 GO TO 905 65 IG = 1 IP = KP - 1 C :::::::::: JACOBI STEP MODIFIED :::::::::: 70 DO 90 K = IG, KP CALL OPB(N,X(1,K),W(1,1)) C DO 85 J = 1, N 85 X(J,K) = W(J,1) C 90 CONTINUE C L1 = IG II = 2 GO TO 905 100 IF (KS .GT. 0) GO TO 130 C :::::::::: MEASURES AGAINST UNHAPPY CHOICE OF C INITIAL VECTORS :::::::::: DO 120 K = 1, KP IF (B(K,K) .NE. 0.0D0) GO TO 120 C DO 110 J = 1, N KZ = MOD(3125*KZ,65536) X(J,K) = DFLOAT(KZ-32768) 110 CONTINUE C KS = 1 120 CONTINUE C IF (KS .NE. 1) GO TO 130 L1 = 1 II = 3 GO TO 905 C 130 DO 150 K = IG, KP C DO 145 L = K, KP T = 0.0D0 C DO 140 I = L, KP 140 T = T + B(I,K) * B(I,L) C :::::::::: NEGATE MATRIX TO REVERSE EIGENVALUE ORDERING :::::::::: B(L,K) = -T 145 CONTINUE C 150 CONTINUE C :::::::::: TRED2, TQL2 ARE MEMBERS OF EISPACK PACKAGE :::::::::: J = KP - KG CALL TRED2(KP,J,B(IG,IG),D(IG),U,B(IG,IG)) CALL TQL2(KP,J,D(IG),U,B(IG,IG),II) C DO 190 K = IG, KP 190 D(K) = DSQRT(DMAX1(-D(K),0.0D0)) C C DO 200 K = IG, KP C DO 200 J = 1, N C 200 W(J,K) = 0.0D0 C C DO 205 K = IG, KP C DO 205 L = IG, KP C DO 205 J = 1, N C 205 W(J,K) = W(J,K) + X(J,L) * B(L,K) C CALL DGEMM('N','N',N,KP-IG+1,KP-IG+1,1.0D0,X(1,IG),NM, * B(IG,IG),KP,0.0D0,W(1,IG),NM) C C DO 210 K = IG, KP DO 210 J = 1, N 210 X(J,K) = W(J,K) C KS = KS + 1 XKS = DFLOAT(KS) IF (D(KP) .GT. E) E = D(KP) C :::::::::: RANDOMIZATION :::::::::: IF (KZ1 .GE. 3) GO TO 225 C DO 220 J = 1, N KZ = MOD(3125*KZ,65536) X(J,KP) = DFLOAT(KZ-32768) 220 CONTINUE C L1 = KP II = 4 GO TO 905 C :::::::::: COMPUTE CONTROL QUANTITIES CX :::::::::: 225 DO 270 K = IG, IP T = (D(K) - E) * (D(K) + E) IF (T .GT. 0.0D0) GO TO 240 CX(K) = 0.0D0 GO TO 270 240 IF (E .NE. 0.0D0) GO TO 260 CX(K) = 1.0D3 + DLOG(D(K)) GO TO 270 260 CX(K) = DLOG((D(K)+DSQRT(T))/E) 270 CONTINUE C C :::::::::: ACCEPTANCE TEST FOR EIGENVALUES INCLUDING ADJUSTMENT C OF KEM AND KH SUCH THAT D(KEM) .GT. E, D(KH) .GT. E, C AND D(KEM) DOES NOT OSCILLATE STRONGLY :::::::::: DO 300 K = IG, KEM IF (D(K) .LE. E .OR. X (KZ1 .GT. 1 .AND. D(K) .LE. 9.99D-1 * RQ(K))) GO TO 320 300 CONTINUE C GO TO 340 320 KEM = K - 1 340 IF (KEM .EQ. 0) GO TO 900 350 K = KH + 1 IF (D(K) .EQ. 0.0D0 .OR. D(K) .GT. EE2 * RQ(K)) GO TO 360 KH = K GO TO 350 360 IF (D(K) .LE. E) KH = K - 1 K = K - 1 IF (K .GT. KEM) GO TO 360 C :::::::::: ACCEPTANCE TEST FOR EIGENVECTORS :::::::::: L = KG E2 = 0.0D0 C DO 560 K = IG, IP IF (K .NE. L + 1) GO TO 430 C :::::::::: CHECK FOR NESTED EIGENVALUES :::::::::: L = K L1 = K IF (K .EQ. IP) GO TO 410 IK = K + 1 S(1) = 5.0D-1 / XKS T = 1.0D0 / DFLOAT(KS*M) C DO 400 J = IK, IP IF (CX(J) * (CX(J) + S(1)) + T .LE. CX(J-1) * CX(J-1)) X GO TO 410 L = J 400 CONTINUE C 410 IF (L .LE. KH) GO TO 430 L = L1 - 1 GO TO 570 430 CALL OPB(N,X(1,K),W(1,1)) S(1) = 0.0D0 C DO 480 J = 1, L IF (DABS(D(J)-D(K)) .GE. 1.0D-2 * D(K)) GO TO 480 T = 0.0D0 C DO 460 I = 1, N 460 T = T + W(I,1) * X(I,J) C DO 470 I = 1, N 470 W(I,1) = W(I,1) - T * X(I,J) C S(1) = S(1) + T * T 480 CONTINUE C T = 0.0D0 C DO 490 I = 1, N 490 T = T + W(I,1) * W(I,1) C IF (S(1) .NE. 0.0D0) GO TO 510 T = 1.0D0 GO TO 520 510 T = DSQRT(T/(S(1)+T)) 520 IF (T .GT. E2) E2 = T IF (K .NE. L) GO TO 560 C :::::::::: TEST FOR ACCEPTANCE OF GROUP OF EIGENVECTORS :::::::::: IF (L .GE. KEM .AND. X D(KEM) * F(KEM) .LT. EPS * (D(KEM) - E)) KG = KEM IF (E2 .GE. F(L)) GO TO 555 C DO 550 J = L1, L 550 F(J) = E2 C 555 IF (L .LE. KEM .AND. D(L) * F(L) .LT. EPS * (D(L) - E)) KG = L IG = KG + 1 560 CONTINUE C :::::::::: ADJUST M :::::::::: 570 IF (E .GT. 4.0D-2 * D(1)) GO TO 590 M = 1 K = 1 GO TO 600 590 E2 = 2.0D0 / E E1 = 5.1D-1 * E2 K = 2 * MAX0(IDINT(4.0D0/CX(1)),1) IF (M .GT. K) M = K C :::::::::: REDUCE KEM IF CONVERGENCE WOULD BE TOO SLOW :::::::::: 600 XKM = DFLOAT(KM) IF (F(KEM) .EQ. 0.0D0 .OR. XKS .GE. 9.0D-1 * XKM) GO TO 660 XK = DFLOAT(K) S(1) = XK * CX(KEM) IF (S(1) .GE. 5.0D-2) GO TO 640 T = 5.0D-1 * S(1) * CX(KEM) GO TO 650 640 T = CX(KEM) + DLOG(5.0D-1*(1.0D0+DEXP(-2.0D0*S(1)))) / XK 650 S(1) = DLOG(D(KEM)*F(KEM)/(EPS*(D(KEM)-E))) / T IF ((XKM - XKS) * XKM .LT. S(1) * XKS) KEM = KEM - 1 660 CALL INF(KS,KG,KH,F,M) C IF (KG .GE. KEM .OR. KS .GE. KM) GO TO 900 C DO 670 K = IG, IP 670 RQ(K) = D(K) C 680 IF (KS + M .LE. KM) GO TO 700 KZ2 = -1 IF (M .GT. 1) M = 2 * ((KM - KS + 1) / 2) 700 M1 = M C :::::::::: SHORTCUT LAST INTERMEDIATE BLOCK IF ALL ERROR C QUANTITIES F ARE SUFFICIENTLY SMALL :::::::::: IF (L .LT. KEM) GO TO 750 S(1) = D(KEM) * F(KEM) / (EPS * (D(KEM) - E)) T = S(1) * S(1) - 1.0D0 IF (T .LE. 0.0D0) GO TO 70 S(1) = DLOG(S(1)+DSQRT(T)) / (CX(KEM) - CX(KH+1)) M1 = 2 * IDINT(5.0D-1*S(1)+1.01D0) IF (M1 .LE. M) GO TO 740 M1 = M GO TO 750 740 KZ2 = -1 750 XM1 = DFLOAT(M1) C :::::::::: CHEBYSHEV ITERATION :::::::::: IF (M .NE. 1) GO TO 790 C DO 780 K = IG, KP CALL OPB(N,X(1,K),W(1,1)) C DO 775 I = 1, N 775 X(I,K) = W(I,1) C 780 CONTINUE C GO TO 860 C :::::::::: DEGREE .NE. ONE :::::::::: C 790 DO 850 K = IG, KP CALL OPB(N,X(1,K),W(1,1)) C DO 810 I = 1, N 810 U(I) = E1 * W(I,1) C CALL OPB(N,U(1),W(1,1)) C DO 820 I = 1, N 820 X(I,K) = E2 * W(I,1) - X(I,K) C IF (M1 .LT. 4) GO TO 850 C DO 840 J = 4, M1, 2 CALL OPB(N,X(1,K),W(1,1)) C DO 830 I = 1, N 830 U(I) = E2 * W(I,1) - U(I) C CALL OPB(N,U,W(1,1)) C DO 835 I = 1, N 835 X(I,K) = E2 * W(I,1) - X(I,K) C 840 CONTINUE C 850 CONTINUE C 860 L1 = IG II = 5 GO TO 905 C :::::::::: DISCOUNTING THE ERROR QUANTITIES F :::::::::: 870 IF (KG .EQ. KH) GO TO 895 IF (M .NE. 1) GO TO 885 C DO 880 K = IG, KH 880 F(K) = F(K) * (D(KH+1) / D(K)) C GO TO 895 885 T = DEXP(-XM1*CX(KH+1)) C DO 890 K = IG, KH S(1) = DEXP(-XM1*(CX(K)-CX(KH+1))) F(K) = S(1) * F(K) * (1.0D0 + T * T) / X (1.0D0 + (S(1) * T) ** 2) 890 CONTINUE C 895 KS = KS + M1 KZ2 = KZ2 - M1 C :::::::::: POSSIBLE REPETITION OF INTERMEDIATE STEPS :::::::::: IF (KZ2 .GE. 0) GO TO 680 KZ1 = KZ1 + 1 KZ2 = 2 * KZ1 M = 2 * M GO TO 70 900 KEM = KG C :::::::::: SOLVE EIGENVALUE PROBLEM OF PROJECTION C OF MATRIX A (MATRIX B FOR SVD PROBLEM) :::::::::: L1 = 1 II = 6 C :::::::::: IN-LINE PROCEDURE FOR EXTENDING ORTHONORMALIZATION C TO ALL KP COLUMNS OF X :::::::::: 905 JP = KP IF (II .EQ. 6) JP = KP - 1 DO 915 I = 1, JP DO 906 K = L1, JP 906 B(K,I) = 0.0D0 C IK = L1 IF (I .LT. L1) GO TO 909 C IK = I + 1 C C DO 907 J = 1, N C 907 B(I,I) = B(I,I) + X(J,I) * X(J,I) C B(I,I) = DSQRT(B(I,I)) C B(I,I)=DSQRT(DDOT(N,X(1,I),1,X(1,I),1)) C T = 0.0D0 IF (B(I,I) .NE. 0.0D0) T = 1.0D0 / B(I,I) C C DO 908 J = 1, N C 908 X(J,I) = T * X(J,I) C CALL DSCAL(N,T,X(1,I),1) C 909 CONTINUE C ILEN=JP-IK+1 C C DO 912 K = IK, JP C DO 910 J = 1, N C 910 B(K,I) = B(K,I) + X(J,I) * X(J,K) C DO 911 J = 1, N C 911 X(J,K) = X(J,K) - B(K,I) * X(J,I) C 912 CONTINUE C CALL DGEMV('T',N,ILEN,1.0D0,X(1,IK),NM, * X(1,I),1,0.0D0,B(IK,I),1) DO 912 K=IK,JP CALL DAXPY(N,-B(K,I),X(1,I),1,X(1,K),1) 912 CONTINUE C 915 CONTINUE C GO TO (65,100,70,225,870,920), II C 920 DO 921 K = 1, IP DO 921 I = K, IP 921 B(I,K) = 0.0D0 C DO 930 K = 1, IP CALL OPB(N,X(1,K),W(1,1)) C DO 925 I = 1, K DO 922 L = 1, N 922 B(K,I) = B(K,I) - X(L,I) * W(L,1) C 925 CONTINUE C 930 CONTINUE C CALL TRED2(KP,IP,B,D,U,B) CALL TQL2(KP,IP,D,U,B,II) C :::::::::: REORDERING OF EIGENVALUES AND EIGENVECTORS ACCORDING C TO THE MAGNITUDES OF THE FORMER :::::::::: DO 934 I = 1, IP IF (I .EQ. IP) GO TO 933 K = I T = D(I) II = I + 1 C DO 931 J = II, IP IF (DABS(D(J)) .LE. DABS(T)) GO TO 931 K = J T = D(J) 931 CONTINUE C IF (K .EQ. I) GO TO 933 D(K) = D(I) D(I) = T C DO 932 J = 1, IP S(J) = B(J,I) B(J,I) = B(J,K) B(J,K) = S(J) 932 CONTINUE C 933 D(I) = -D(I) 934 CONTINUE C C USE LIBRARY KERNEL FOR MATRIX MULT. C C DO 940 I = 1, IP C DO 940 J = 1, N C 940 W(J,I) = 0.0D0 C C DO 945 I = 1, IP C DO 945 K = 1, IP C DO 945 J = 1, N C 945 W(J,I) = W(J,I) + X(J,K) * B(K,I) C CALL DGEMM('N','N',N,IP,IP,1.0D0,X,NM,B,KP,0.0D0,W,NM) C DO 950 I = 1, IP DO 950 J = 1, N 950 X(J,I) = W(J,I) C D(KP) = E RETURN C :::::::::: LAST CARD OF RITZIT :::::::::: 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 REAL*8 FUNCTION PYTHAG(A,B) REAL*8 A,B C C FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW C REAL*8 P,R,S,T,U C 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 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 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 SUBROUTINE OPB(N, X, Y) C C.....MULTIPLICATION OF A SHIFTED 2-CYCLIC MATRIX B BY A VECTOR X , WHERE C C B = [ D A ] C [ A' D ] , WHERE A IS NROW BY NCOL (NROW >> NCOL), C AND D = ALPHA*I, ALPHA AN UPPER BOUND C FOR THE LARGEST SINGULAR VALUE OF A. C C HENCE, B IS OF ORDER N:= NROW+NCOL (Y STORES PRODUCT VECTOR) C PARAMETER(NMAX=500,NZMAX=2000) REAL*8 VALUE(NZMAX), ALPHA INTEGER POINTR(NMAX), ROWIND(NZMAX) COMMON /MATRIXA/ NCOL,NROW,VALUE,POINTR,ROWIND,ALPHA C COMMON /COUNT/ MXVCOUNT INTEGER N,MXVCOUNT REAL*8 X(1), Y(1) C C--------------------- C MXVCOUNT=MXVCOUNT+2 DO 55 I=1,N Y(I)=ALPHA*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 C SUBROUTINE OPC(N, X, Y) C C.....MULTIPLICATION OF THE 2-CYCLIC MATRIX B BY A VECTOR X , WHERE C C B = [ 0 A ] C [ A' 0 ] , WHERE A IS NROW BY NCOL (NROW >> NCOL). C C C HENCE, B IS OF ORDER N:= NROW+NCOL (Y STORES PRODUCT VECTOR) C PARAMETER(NMAX=500,NZMAX=2000) REAL*8 VALUE(NZMAX), ALPHA INTEGER POINTR(NMAX), ROWIND(NZMAX) COMMON /MATRIXA/ NCOL,NROW,VALUE,POINTR,ROWIND,ALPHA C COMMON /COUNT/ MXVCOUNT INTEGER N,MXVCOUNT REAL*8 X(1), Y(1) C C--------------------- C MXVCOUNT=MXVCOUNT+2 DO 55 I=1,N Y(I)=0.0D0 55 CONTINUE C DO 15 I=1,NCOL C DO 10 J=POINTR(I),POINTR(I+1)-1 Y(ROWIND(J))=Y(ROWIND(J))+VALUE(J)*X(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