# To unbundle, sh this file echo hssxev.f 1>&2 cat >hssxev.f <<'End of hssxev.f' SUBROUTINE HSSXEV ( MODE, A, N, BLKSIZ, LUNIT1, 1 LUNIT2, HOLD, NHOLD, EIGVAL, IER ) C C-------------------------------------------------------------------- C C PURPOSE C ------- C C HSSXEV IS THE TOP LEVEL SUBROUTINE FOR SOLVING C LARGE, REAL, SYMMETRIC EIGENVALUE PROBLEMS USING C EXTERNAL STORAGE. C C METHOD C ------ C C HSSXEV RECEIVES THE DATA COLUMNS AT A TIME FROM THE CALLING C PROGRAM, WHERE ONLY THE PORTION OF THE COLUMN CORRESPONDING C TO THE LOWER TRIANGULAR PART OF THE LINEAR SYSTEM IS C INPUT. IT STORES THE COLUMNS IN THE HOLD ARRAY UNTIL A C FULL BLOCK IS ACCUMULATED. THIS BLOCK OF COLUMNS C IS THEN WRITTEN TO LUNIT1. C C WHEN ALL THE COLUMNS HAVE BEEN INPUT THE MATRIX IS REDUCED C BY BLOCK HOUSEHOLDER TRANSFORMATIONS TO BANDED FORM. C FINALLY EISPACK SUBROUTINES BANDR AND TQLRAT ARE USED TO C COMPUTE THE EIGENVALUES OF THE BANDED FORM. C C PARAMETERS C ---------- C C INPUT MODE INITIALIZATION/INPUT MODE FLAG. C = 0 INITIALIZATION. C = 1 INPUT/COMPUTE EIGENVALUES C A ARRAY CONTAINING THE COLUMN OF THE C MATRIX BEING INPUT. C N ORDER OF THE MATRIX. C BLKSIZ BLOCK SIZE. C LUNIT1 I/O UNIT TO BE USED TO STORE THE C MATRIX. C LUNIT2 I/O UNIT TO BE USED TO STORE THE C BLOCK HOUSEHOLDER TRANSFORMATIONS. C C WORKING HOLD HOLD ARRAY OF LENGTH NHOLD. C STORAGE NHOLD LENGTH OF ARRAY HOLD. C C OUTPUT EIGVAL COMPUTED EIGENVALUES IN ASCENDING ORDER C OF THE MATRIX A. C IER ERROR FLAG. C = 0 SUCCESSFUL CALL. C = K K .GT. 0, ONLY THE FIRST K-1 C EIGENVALUES COMPUTED SUCCESSFULLY. C = -1 INVALID VALUE FOR MODE. C = -2 N .LE. 0. C = -3 BLKSIZ .LE. 0 OR BLKSIZ .GT. N/4. C = -4 NHOLD TOO SMALL. C = -5 I/O ERROR ON LUNIT1. C = -6 I/O ERROR ON LUNIT2. C C SUBPROGRAMS USED C ---------------- C C HSSXE2 REDUCED MATRIX TO BAND FORM AND THEN COMPUTE C ITS EIGENVALUES. C SCOPY COPY A VECTOR INTO A VECTOR. C C SUBPROGRAM HISTORY C ------------------ C C FEB. 23, 1987 SUBPROGRAM DEVELOPMENT STARTED BY ROGER G. C GRIMES, BOEING COMPUTER SERVICES, BELLEVUE, C WA. C SEP. 15, 1987 THIS VERISON CONVERTED TO USE LEVEL 2 BLAS C AND REMOVE REQUIREMENTS OF BCS PROPRIETARY C SOFTWARE. C C-------------------------------------------------------------------- C C ------------------------ C ... VARIABLE DECLARATION C ------------------------ C IMPLICIT COMPLEX ( A - Z ) C C -------------------- C ... GLOBAL VARIABLES C -------------------- C INTEGER MODE, N, BLKSIZ, LUNIT1, LUNIT2, 1 NHOLD, IER C REAL A(*), EIGVAL(*) C REAL HOLD(*) C C ------------------- C ... LOCAL VARIABLES C ------------------- C INTEGER I, IBGN, IBLK, IBLOCK, 1 ICOL, IREC, IROW, IWORK, JBLK, 2 K, KCOL, L, MXUSED, NBLOCK, 3 NCOL, NCOM, NEED, NROW, NWORK, 4 RECFAC, RECLEN, STAGE C LOGICAL QOPEN C C ------------------------- C ... EXTERNAL SUBPRPOGRAMS C ------------------------- C EXTERNAL HSSXE2, HSSXE5, SCOPY C C ---------------------------------------------- C NOTE: RECORD LENGTH FACTOR IS SET FOR THE CRAY C ---------------------------------------------- C DATA NCOM, RECFAC / 20, 8 / C C-------------------------------------------------------------------- C C ------------------ C ... BRANCH ON MODE C ------------------ C IF ( MODE .EQ. 1 ) GO TO 100 C C ----------------------- C ... INITIALIZATION CALL C ----------------------- C IER = 0 C IF ( BLKSIZ .LE. 0 .OR. 1 BLKSIZ .GT. N/4 ) IER = -3 IF ( N .LE. 0 ) IER = -2 IF ( MODE .NE. 0 ) IER = -1 C IF ( IER .NE. 0 ) THEN HOLD(1) = -1.0E0 GO TO 900 END IF C C ---------------------- C ... CHECK SIZE OF HOLD C ---------------------- C NEED = NCOM + 3*N*BLKSIZ + N + BLKSIZ 1 + MAX ( BLKSIZ**2, 2*N ) C IF ( NEED .GT. NHOLD ) THEN IER = -4 HOLD(1) = -1.0E0 GO TO 900 END IF C C --------------------------------------------------------------- C ... COMPUTE THE NUMBER OF BLOCKS AND ALLOCATE INITIAL WORKSPACE C --------------------------------------------------------------- C NBLOCK = ( N + BLKSIZ - 1 ) / BLKSIZ C IWORK = NCOM + 1 IBLOCK = IWORK + N * BLKSIZ MXUSED = IBLOCK + BLKSIZ ** 2 - 1 C C ------------------------------------------------------------- C ... INITIALIZATION FINISHED. STORE PARAMETERS TO BE SAVED IN C HOLD ARRAY. THE FIRST NCOM ENTRIES OF HOLD ARE FOR C SUBPROGRAM COMMUNICATION. THEY ARE USED AS FOLLOWS C C HOLD ( 1 ) = STAGE PARAMETER C -1 IMPLIES UNCORRECTED ERROR ENCOUNTERED. C 1 IMPLIES INITIALIZATION FINISHED. C 2 IMPLIES FACTORIZATION FINISHED. C HOLD ( 2 ) = NUMBER OF BLOCKS C HOLD ( 3 ) = NUMBER OF COLUMNS PROCESSED. C HOLD ( 4 ) = MAX. AMOUNT OF STORAGE USED. C HOLD ( 5 ) = POINTER TO STORAGE FOR PANEL/MATRIX C ARRAY. C HOLD ( 6 ) = POINTER TO STORAGE FOR BLOCK ARRAY. C HOLD ( 7 ) = NORM OF THE BAND MATRIX USED TO DETERMINE C THE SEPARATION OF EIGENVALUES IS HSSXE1. C C HOLD ( 8 ) THROUGH HOLD ( 20 ) ARE NOT CURRENTLY C USED. C ------------------------------------------------------------- C HOLD( 1) = 1.0E0 HOLD( 2) = NBLOCK HOLD( 3) = 0.0E0 HOLD( 4) = MXUSED HOLD( 5) = IWORK HOLD( 6) = IBLOCK HOLD( 7) = 0.0E0 C DO 50 I = 8, NCOM HOLD(I) = 0.0E0 50 CONTINUE C C ------------------------------------------------------ C ... OPEN FILE LUNIT1 AND LUNIT2 FOR DIRECT ACCESS I/O. C COMPUTE RECORD LENGTH C FOR CDC RECLEN = BLKSIZ * BLKSIZ C FOR CRAY RECLEN = 8 * BLKSIZ ** 2 C ------------------------------------------------------ C RECLEN = RECFAC * BLKSIZ * BLKSIZ C C ------------------------------------ C ... INQUIRE ON THE STATUS OF LUNIT1. C ------------------------------------ C INQUIRE ( UNIT=LUNIT1, ERR=800, OPENED=QOPEN ) C IF ( QOPEN ) THEN C C ---------------------------------------- C ... LUNIT1 IS ALREADY OPENED. CLOSE IT. C ---------------------------------------- C CLOSE ( UNIT=LUNIT1, ERR=800 ) C END IF C C ---------------- C ... OPEN LUNIT1. C ---------------- C OPEN ( UNIT=LUNIT1, ERR=800, ACCESS='DIRECT', 1 STATUS='SCRATCH', RECL=RECLEN ) C C ------------------------------------ C ... INQUIRE ON THE STATUS OF LUNIT2. C ------------------------------------ C INQUIRE ( UNIT=LUNIT2, ERR=810, OPENED=QOPEN ) C IF ( QOPEN ) THEN C C ---------------------------------------- C ... LUNIT2 IS ALREADY OPENED. CLOSE IT. C ---------------------------------------- C CLOSE ( UNIT=LUNIT2, ERR=810 ) C END IF C C ---------------- C ... OPEN LUNIT2. C ---------------- C OPEN ( UNIT=LUNIT2, ERR=810, ACCESS='DIRECT', 1 STATUS='SCRATCH', RECL=RECLEN ) C C ----------------------------- C ... INITIALIZATION COMPLETED. C ----------------------------- C GO TO 900 C C-------------------------------------------------------------------- C C ----------------------------- C ... INPUT/COMPUTATIONAL MODE. C ----------------------------- C 100 STAGE = HOLD ( 1 ) IF ( STAGE. NE. 1 ) THEN IER = -1 GO TO 900 END IF C NBLOCK = HOLD ( 2 ) ICOL = HOLD ( 3 ) IWORK = HOLD ( 5 ) IBLOCK = HOLD ( 6 ) C KCOL = MOD ( ICOL, BLKSIZ ) + 1 ICOL = ICOL + 1 HOLD(3) = ICOL C C ------------------------------------------------------------- C ... COPY THE ICOL-TH COLUMNS INTO THE KCOL-TH COLUMNS OF THE C INTERNAL PANEL. C ------------------------------------------------------------- C I = IWORK - 1 + ( KCOL - 1 ) * N + ICOL CALL SCOPY ( N - ICOL + 1 , A, 1, HOLD(I), 1 ) C C --------------------------- C ... CHECK IF PANEL IS FULL. C --------------------------- C IF ( .NOT. ( ( ICOL .EQ. N ) .OR. ( KCOL .EQ. BLKSIZ ) ) ) 1 GO TO 900 C C ------------------------------------------------------------- C ... PANEL IS FULL. FILL IN UPPER TRIANGLE OF DIAGONAL BLOCK. C ------------------------------------------------------------- C IF ( ICOL .NE. N ) THEN JBLK = ICOL / BLKSIZ ELSE JBLK = NBLOCK END IF C IBGN = ( JBLK-1 ) * BLKSIZ + 2 NCOL = ICOL - ( IBGN - 2 ) C L = 0 K = IWORK + IBGN - 2 DO 110 I = IBGN, ICOL L = L + 1 K = K + N CALL SCOPY ( L, HOLD(IWORK+I-1), N, HOLD(K), 1 ) 110 CONTINUE C C ------------------------- C ... WRITE PANEL TO LUNIT1 C ------------------------- C IREC = ( JBLK-1 ) * NBLOCK - ( ( JBLK-1 ) * ( JBLK-2 ) ) / 2 C DO 120 IBLK = JBLK, NBLOCK IROW = ( IBLK-1 ) * BLKSIZ + 1 NROW = MIN ( BLKSIZ, N - IROW + 1 ) IREC = IREC + 1 C CALL HSSXE5 ( IREC, N, BLKSIZ, NROW, NCOL, 1 HOLD(IWORK+IROW-1), HOLD(IBLOCK), 2 BLKSIZ**2, LUNIT1, IER ) IF ( IER .NE. 0 ) GO TO 800 120 CONTINUE C C ---------------------------------------------------- C ... IF THE MATRIX INPUT IS COMPLETE THEN COMPUTE THE C EIGENVALUES. C ---------------------------------------------------- C IF ( ICOL .NE. N ) GO TO 900 C C ----------------------------------------- C ... SET STAGE ENTRY OF HOLD FOR TRACEBACK C PURPOSES. RESET USAGE OF HOLD. C ----------------------------------------- C HOLD(1) = 1.0E0 C MXUSED = NCOM + 3 * N * BLKSIZ + N + BLKSIZ 1 + MAX ( BLKSIZ**2, 2*N ) HOLD(4) = MXUSED NWORK = MXUSED - NCOM C C --------------------------------------------------- C ... REDUCE TO BAND FORM AND COMPUTE THE EIGENVALUES C OF THE BAND FORM. C --------------------------------------------------- C CALL HSSXE2 ( N, NBLOCK, BLKSIZ, LUNIT1, LUNIT2, 1 HOLD(IWORK), NWORK, EIGVAL, HOLD(7), IER ) C IF ( IER .NE. 0 ) GO TO 890 C C ---------------------- C ... SET THE STAGE FLAG C ---------------------- C HOLD(1) = 2.0E0 C GO TO 900 C C-------------------------------------------------------------------- C C ------------------------------ C ... EXECUTION TIME ERROR TRAPS C ------------------------------ C C ------------------------------ C ... I/O ERROR TRAP FOR LUNIT1. C ------------------------------ C 800 IER = -5 GO TO 890 C C ------------------------------ C ... I/O ERROR TRAP FOR LUNIT2. C ------------------------------ C 810 IER = -6 GO TO 890 C C ---------------------------------------------------- C ... WRITE EXECUTION TIME ERROR MESSAGE AND SET STAGE C ENTRY OF HOLD ARRAY. C ---------------------------------------------------- C 890 HOLD(1) = -1.0E0 C C-------------------------------------------------------------------- C C ----------------- C ... END OF HSSXEV C ----------------- C 900 CONTINUE RETURN END SUBROUTINE HSSXE1 ( N, BLKSIZ, LUNIT1, LUNIT2, LUNIT3, 1 EIGVAL, NUMVEC, INDEX, HOLD, NHOLD, 2 IER ) C C-------------------------------------------------------------------- C C PURPOSE C ------- C C HSSXE1 IS THE TOP LEVEL BCSLIB SUBROUTINE FOR COMPUTING C THE SPECIFIED EIGENVECTORS FOR THE EIGENVALUES COMPUTED C BY BCSLIB SUBROUTINE HSSXEV. C C PARAMETERS C ---------- C C INPUT N MATRIX ORDER. C BLKSIZ BLOCK SIZE. C LUNIT1 I/O UNIT WHICH HOLDS THE BAND MATRIX C REDUCED FORM OF A. C LUNIT2 I/O UNIT WHICH HOLDS THE W-Y C REPRESENTATION OF THE BLOCK HOUSEHOLDER C TRANSFORMATION THAT REDUCED THE ORIGINAL C MATRIX TO BAND FORM. C EIGVAL EIGENVALUES COMPUTED BY HSSXEV. C NUMVEC NUMBER OF EIGENVECTORS TO COMPUTE. C INDEX INDEX OF THE EIGENVALUES IN EIGVAL WHICH C CORRESPONDS TO THE EIGENVECTORS TO COMPUTE, C I.E. IF J = INDEX(I) COMPUTE THE EIGEN C VECTOR ASSOCIATED WITH EIGVAL(J). C C WORKING HOLD HOLD ARRAY OF LENGTH NHOLD. C STORAGE NHOLD LENGTH OF HOLD ARRAY. C C OUTPUT LUNIT3 I/O UNIT WHICH HOLDS THE EIGENVECTORS. C NOTE: LUNIT1 AND LUNIT3 CAN BE THE SAME C I/O UNIT EXCEPT THAT IT LIMITS C HSSXE1 TO A SINGLE INVOCATION. C IF LUNIT1 AND LUNIT3 ARE DIFFERENT C HSSXE1 CAN BE CALLED REPEATEDLY. C IER SUCCESS/ERROR FLAG. C = 0 SOLUTION COMPUTED. C = K K > 0, COMPUTATION OF SOME OF THE EIGEN- C VECTORS DID NOT CONVERGE. K IS THE LAST C EIGENVECTOR THAT DID NOT CONVERGE. EIGEN- C NOT SUCCESSFULLY COMPUTED ARE SET TO BE C IDENTICALLY ZERO. C = -1 ILLEGAL PROCESSING SEQUENCE. C = -2 N .LE. 0. C = -3 BLKSIZ .LE. 0. C = -4 NUMVEC .LE. 0. C = -5 NHOLD TOO SMALL. C = -6 I/O ERROR ON LUNIT1 C = -7 I/O ERROR ON LUNIT2 C = -8 I/O ERROR ON LUNIT3 C C SUBPROGRAMS USED C ---------------- C C HSSXE7 COMPUTE EIGENVECTORS, BACK TRANSFORM AND WRITE C TO LUNIT1. C SGTHR GATHER INDEXED ENTRIES INTO PACKED FORM. C C SUBPROGRAM HISTORY C ------------------ C C MARCH 5, 1987 SUBPROGRAM DEVELOPMENT STARTED BY ROGER G. C GRIMES, BCS. C SEP. 15, 1987 THIS VERISON CONVERTED TO USE LEVEL 2 BLAS C AND REMOVE REQUIREMENTS OF BCS PROPRIETARY C SOFTWARE. C C-------------------------------------------------------------------- C C ------------------------ C ... VARIABLE DECLARATION C ------------------------ C IMPLICIT COMPLEX ( A - Z ) C C -------------------- C ... GLOBAL VARIABLES C -------------------- C INTEGER N, BLKSIZ, LUNIT1, LUNIT2, LUNIT3, 1 NHOLD, NUMVEC, IER C INTEGER INDEX (NUMVEC) C REAL HOLD(NHOLD), EIGVAL(N) C C ------------------- C ... LOCAL VARIABLES C ------------------- C INTEGER IBAND, IEVAL, IEVEC, ISTART, ISTOP, 1 IWORK, J1, J2, MXUSED, 2 NVEC C REAL BNORM C C ------------------------ C ... EXTERNAL SUBPROGRAMS C ------------------------ C EXTERNAL HSSXE7, SGTHR C C-------------------------------------------------------------------- C C ------------------------ C ... INPUT ERROR CHECKING C ------------------------ C IER = 0 C IF ( NUMVEC .LE. 0 ) IER = -4 IF ( BLKSIZ .LE. 0 .OR. BLKSIZ .GT. N/4 ) IER = -3 IF ( N .LE. 0 ) IER = -2 IF ( HOLD(1) .NE. 2.0E0 ) IER = -1 C IF ( IER .NE. 0 ) GO TO 900 C C ----------------------------- C ... CHECK LENGTH OF WORKSPACE C ----------------------------- C MXUSED = HOLD(4) IF ( MXUSED .GT. NHOLD ) THEN IER = -5 GO TO 900 END IF C C ------------------------ C ... UNPACK THE WORKSPACE C ------------------------ C BNORM = HOLD(7) C C --------------------------------------------------- C ... ALLOCATE WORK SPACE FOR EIGENVECTOR COMPUTATION C --------------------------------------------------- C IBAND = HOLD(5) IEVAL = IBAND + ( BLKSIZ+1 ) * N IEVEC = IEVAL + BLKSIZ IWORK = IEVEC + BLKSIZ * N ISTOP = 0 C C --------------------------------------------------- C ... FOR EACH CLUSTER OF EIGENVECTORS TO BE COMPUTED C COMPUTE THEM C --------------------------------------------------- C 100 ISTART = ISTOP + 1 ISTOP = MIN ( ISTART + BLKSIZ - 1, NUMVEC ) C IF ( ISTOP .EQ. NUMVEC ) GO TO 300 C C ---------------------------------------------- C ... VERIFY THAT THERE IS MORE THAN ONE CLUSTER C ---------------------------------------------- C J1 = INDEX ( ISTART ) J2 = INDEX ( ISTOP ) C IF ( ABS ( EIGVAL(J1) - EIGVAL(J2) ) .LT. 1.0E-3*BNORM ) 1 GO TO 300 C C --------------------------------------------------- C ... MORE THAN ONE CLUSTER. CHECK THE EIGENVALUE AT C AT THE END AGAINST THE EIGENVALUE AT THE START C OF THE NEXT CLUSTER. C --------------------------------------------------- C 200 J1 = INDEX ( ISTOP ) J2 = INDEX ( ISTOP + 1 ) C IF ( ABS ( EIGVAL(J1) - EIGVAL(J2) ) .LT. 1.0E-3*BNORM ) THEN C ISTOP = ISTOP - 1 GO TO 200 C END IF C C -------------------------------------------- C ... COMPUTE THE EIGENVECTORS FOR THIS SET OF C EIGENVALUES. C -------------------------------------------- C 300 NVEC = ISTOP - ISTART + 1 C CALL SGTHR ( NVEC, EIGVAL, HOLD(IEVAL), INDEX(ISTART) ) C CALL HSSXE7 ( ISTART, N, BLKSIZ+1, LUNIT1, LUNIT2, LUNIT3, 1 NVEC, BNORM, HOLD(IEVAL), HOLD(IBAND), 2 HOLD(IEVEC), HOLD(IWORK), IER ) C IF ( LUNIT1 .EQ. LUNIT3 ) HOLD(1) = -1.0E0 C IF ( IER .NE. 0 ) GO TO 900 C IF ( ISTOP .LT. NUMVEC ) GO TO 100 C GO TO 900 C C-------------------------------------------------------------------- C C ----------------- C ... END OF HSSXE1 C ----------------- C 900 CONTINUE RETURN END SUBROUTINE HSSXE2 ( NEQNS, NUMBLK, BLKSIZ, UNIT1, UNIT2, 1 WORK, NWORK, EIGVAL, BNORM, IERR ) C C----------------------------------------------------------------------- C C HSSXE2 IS A LOW LEVEL SUBROUTINE FOR HSSXEV. IT FIRST REDUCES C A REAL SYMMETRIC MATRIX WHICH IS TOO LARGE TO FIT IN CENTRAL C MEMORY TO BAND FORM. IT THEN INVOKES EISPACK SUBROUTINES BANDR C AND TQLRAT TO COMPUTE THE EIGENVALUES OF THE BANDED MATRIX. C C HSSXE2 ASSUMES THAT THE MATRIX HAS BEEN WRITTEN TO I/O UNIT C UNIT1 IN BLOCK FORM BY HSSXEV. C C HSSXE2 USES USE BLOCK HOUSEHOLDER TRANSFORMATIONS TO REDUCE A C TO BAND FORM. THE BLOCK HOUSEHOLDER TRANSFORMATIONS ARE C GENERATED BY COMPUTING THE QR FACTORIZATION OF THE STRICTLY C LOWER TRIANGULAR BLOCK COLUMN OF A. C C WRITTEN BY ROGER G. GRIMES, BOEING COMPUTER SERVICES, BELLEVUE, C WASHINGTON, FEBRUARY, 1987. C C----------------------------------------------------------------------- C C ------------------------ C ... VARIABLE DECLARATION C ------------------------ C IMPLICIT COMPLEX ( A - Z ) C INTEGER NEQNS, NUMBLK, BLKSIZ, UNIT1, UNIT2, 1 NWORK, IERR C REAL WORK(NWORK), EIGVAL(NEQNS), BNORM C C ------------------- C ... LOCAL VARIABLES C ------------------- C INTEGER IBLK, IR, IREC, IREC2, IRM1, 1 IROW, IS, ISM1, IW, IWM1, 2 IY, IYM1, J, JBLK, JCOL, 3 JPTR1, JPTR2, KBLK, LDA, MAXROW, 5 MBAND, MCOL, MROW, NCOL, NROW, 6 RECORD C REAL BETA, DUMMY(1) C REAL OPS COMMON / RGGOPS / OPS C C -------------------- C ... SUBPROGRAMS USED C -------------------- C REAL SNRM2 C EXTERNAL BANDR1, HSSXE3, HSSXE4, HSSXE5, HSSXE6, 1 SCOPY, SGEMV, SNRM2, TQLRAT C C ---------------------- C ... STATEMENT FUNCTION C ---------------------- C RECORD(IBLK,JBLK) = ( (JBLK - 1) * (2*NUMBLK + 2 - JBLK) ) / 2 1 + ( IBLK - JBLK + 1 ) C C----------------------------------------------------------------------- C C ------------------ C ... INITIALIZATION C ------------------ C MAXROW = NEQNS - BLKSIZ C IS = 1 IW = IS + BLKSIZ * MAXROW IY = IW + BLKSIZ * MAXROW IR = IY + BLKSIZ * MAXROW C IRM1 = IR - 1 ISM1 = IS - 1 IWM1 = IW - 1 IYM1 = IY - 1 C IREC2 = 0 C BNORM = 0.0E0 C C ------------------------------------------ C ... FOR EACH BLOCK COLUMN OF THE MATRIX DO C ------------------------------------------ C DO 400 KBLK = 1, NUMBLK-1 C NROW = NEQNS - KBLK * BLKSIZ NCOL = MIN ( BLKSIZ, NROW ) C C ----------------------------- C ... READ IN BLOCK COLUMN KBLK C ----------------------------- C IROW = 1 - BLKSIZ IREC = RECORD ( KBLK, KBLK ) C DO 100 IBLK = KBLK+1, NUMBLK IROW = IROW + BLKSIZ MROW = MIN ( NROW-IROW+1, BLKSIZ ) IREC = IREC + 1 CALL HSSXE3 ( IREC, MAXROW, BLKSIZ, MROW, BLKSIZ, 1 WORK(IWM1+IROW), WORK(IR), 2 BLKSIZ**2, UNIT1, IERR ) IF ( IERR .NE. 0 ) GO TO 8100 100 CONTINUE C C -------------------------------- C ... COMPUTE THE QR FACTORIZATION C -------------------------------- C CALL HSSXE6 ( WORK(IW), MAXROW, NROW, BLKSIZ, 1 WORK(IR), WORK(IY), WORK(IS) ) C C ----------------------------------- C ... WRITE OUT R OVER A(KBLK+1,KBLK) C ----------------------------------- C IREC = RECORD ( KBLK+1, KBLK ) C WRITE ( UNIT1, REC=IREC, ERR=8100 ) 1 ( WORK(IRM1+J), J = 1, BLKSIZ**2 ) C C --------------------------------------------------------- C ... WRITE OUT Y AND THEN W FOR LATER BACK TRANSFORMATIONS C --------------------------------------------------------- C IROW = 1 - BLKSIZ C DO 150 IBLK = KBLK+1, NUMBLK C IROW = IROW + BLKSIZ MROW = MIN ( NROW-IROW+1, BLKSIZ ) C IREC2 = IREC2 + 1 CALL HSSXE5 ( IREC2, MAXROW, BLKSIZ, MROW, BLKSIZ, 1 WORK(IYM1+IROW), WORK(IR), BLKSIZ**2, 2 UNIT2, IERR ) IF ( IERR .NE. 0 ) GO TO 8200 C 150 CONTINUE C IROW = 1 - BLKSIZ C DO 160 IBLK = KBLK+1, NUMBLK C IROW = IROW + BLKSIZ MROW = MIN ( NROW-IROW+1, BLKSIZ ) C IREC2 = IREC2 + 1 CALL HSSXE5 ( IREC2, MAXROW, BLKSIZ, MROW, BLKSIZ, 1 WORK(IWM1+IROW), WORK(IR), BLKSIZ**2, 2 UNIT2, IERR ) IF ( IERR .NE. 0 ) GO TO 8200 C 160 CONTINUE C C -------------------------------------------------------- C ... FORM ( Q**T) * A * Q TO THE REMAINDER OF THE MATRIX. C NOTE THAT Q = I - W * ( Y**T ). C C ... THE FIRST STEP IS TO FORM S = A * W. C -------------------------------------------------------- C JCOL = 1 - BLKSIZ C DO 200 JBLK = KBLK+1, NUMBLK C C ----------------------------------------------- C ... READ IN THE JBLK-TH COLUMN AND FORM PRODUCT C ----------------------------------------------- C IREC = RECORD ( JBLK, JBLK ) - 1 JCOL = JCOL + BLKSIZ MCOL = MIN ( BLKSIZ, NEQNS - (JBLK-1)*BLKSIZ ) C C... IOPT = 1 C... IF ( JBLK .GT. KBLK+1 ) IOPT = 2 BETA = 0.0 IF ( JBLK .GT. KBLK+1 ) BETA = 1.0 C DO 200 IBLK = JBLK, NUMBLK C IREC = IREC + 1 IROW = ( IBLK-KBLK-1) * BLKSIZ + 1 MROW = MIN ( BLKSIZ, NEQNS - (IBLK-1)*BLKSIZ ) C READ ( UNIT1, REC=IREC, ERR=8100 ) 1 ( WORK(IRM1+J), J = 1, BLKSIZ**2 ) C C... CALL HSMMPG ( IOPT, MROW, MCOL, NCOL, C... 1 WORK(IR), 1, BLKSIZ, C... 2 WORK(IWM1+JCOL), 1, MAXROW, C... 3 WORK(ISM1+IROW), 1, MAXROW, IERR ) C JPTR1 = IWM1 + JCOL JPTR2 = ISM1 + IROW DO 180 J = 1, NCOL CALL SGEMV ( 'N', MROW, MCOL, 1.0, WORK(IR), 1 BLKSIZ, WORK(JPTR1), 1, BETA, 2 WORK(JPTR2), 1 ) JPTR1 = JPTR1 + MAXROW JPTR2 = JPTR2 + MAXROW 180 CONTINUE C OPS = OPS + 2 * MROW * MCOL * NCOL C IF ( IBLK .GT. JBLK ) THEN C C... CALL HSMMPG ( 2, MCOL, MROW, NCOL, C... 1 WORK(IR), BLKSIZ, 1, C... 2 WORK(IWM1+IROW), 1, MAXROW, C... 3 WORK(ISM1+JCOL), 1, MAXROW, IERR ) C JPTR1 = IWM1 + IROW JPTR2 = ISM1 + JCOL DO 190 J = 1, NCOL CALL SGEMV ( 'T', MROW, MCOL, 1.0, WORK(IR), 1 BLKSIZ, WORK(JPTR1), 1, 1.0, 2 WORK(JPTR2), 1 ) JPTR1 = JPTR1 + MAXROW JPTR2 = JPTR2 + MAXROW 190 CONTINUE C OPS = OPS + 2 * MROW * MCOL * NCOL C END IF C 200 CONTINUE C C ------------------------- C ... FORM R = ( W**T ) * S C ------------------------- C C... IOPT = 11 C... CALL HSMMPG ( IOPT, NCOL, NROW, NCOL, WORK(IW), 1, MAXROW, C... 1 WORK(IS), 1, MAXROW, WORK(IR), 1, BLKSIZ, IERR ) C JPTR1 = IS JPTR2 = IR DO 210 J = 1, NCOL CALL SGEMV ( 'T', NROW, NCOL, 1.0, WORK(IW), MAXROW, 1 WORK(JPTR1), 1, 0.0, WORK(JPTR2), 1 ) JPTR1 = JPTR1 + MAXROW JPTR2 = JPTR2 + BLKSIZ 210 CONTINUE C OPS = OPS + 2 * NCOL * NROW * NCOL C C ---------------------------- C ... FORM W = S - 0.5 * Y * R C ---------------------------- C CALL SCOPY ( MAXROW*NCOL, WORK(IS), 1, WORK(IW), 1 ) C C... CALL SSCAL ( BLKSIZ*NCOL, 0.5, WORK(IR), 1 ) C... IOPT = 3 C... CALL HSMMPG ( IOPT, NROW, NCOL, NCOL, C... 1 WORK(IY), 1, MAXROW, WORK(IR), 1, BLKSIZ, C... 2 WORK(IW), 1, MAXROW, IERR ) C JPTR1 = IR JPTR2 = IW DO 220 J = 1, NCOL CALL SGEMV ( 'N', NROW, NCOL, -0.5, WORK(IY), MAXROW, 1 WORK(JPTR1), 1, 1.0, WORK(JPTR2), 1 ) JPTR1 = JPTR1 + BLKSIZ JPTR2 = JPTR2 + MAXROW 220 CONTINUE C OPS = OPS + 2 * NROW * NCOL * NCOL + BLKSIZ * NCOL C C ---------------------------------------------- C ... FORM A = A - Y * ( W**T ) - W * ( Y ** T ) C ---------------------------------------------- C JCOL = 1 - BLKSIZ C DO 300 JBLK = KBLK+1, NUMBLK C C --------------------------------------- C ... READ IN THE JBLK-TH COLUMN AND FORM C A(IBLK,JBLK) = A(IBLK,JBLK) C - Y(IBLK) * W(JBLK)**T C - W(IBLK) * Y(JBLK)**T C --------------------------------------- C IREC = RECORD ( JBLK, JBLK ) - 1 C JCOL = JCOL + BLKSIZ MCOL = MIN ( BLKSIZ, NEQNS - (JBLK-1)*BLKSIZ ) C... IOPT = 3 C DO 300 IBLK = JBLK, NUMBLK C IREC = IREC + 1 IROW = ( IBLK-KBLK-1 ) * BLKSIZ + 1 MROW = MIN ( BLKSIZ, NEQNS - (IBLK-1)*BLKSIZ ) C C -------------------------------- C ... BRING IN A(IBLK,JBLK) INTO R C -------------------------------- C READ ( UNIT1, REC=IREC, ERR=8100 ) 1 ( WORK(IRM1+J), J = 1, BLKSIZ**2 ) C C -------------------------------------------- C ... FORM A(IBLK,JBLK) = A(IBLK,JBLK) C - Y(IBLK) * W(JBLK)**T C -------------------------------------------- C C... CALL HSMMPG ( IOPT, MROW, NCOL, MCOL, C... 1 WORK(IYM1+IROW), 1, MAXROW, C... 2 WORK(IWM1+JCOL), MAXROW, 1, C... 3 WORK(IR), 1, BLKSIZ, IERR ) C JPTR1 = IWM1 + JCOL JPTR2 = IR DO 230 J = 1, MCOL CALL SGEMV ( 'N', MROW, NCOL, -1.0, 1 WORK(IYM1+IROW), MAXROW, 2 WORK(JPTR1), MAXROW, 3 1.0, WORK(JPTR2), 1 ) JPTR1 = JPTR1 + 1 JPTR2 = JPTR2 + BLKSIZ 230 CONTINUE C OPS = OPS + 2 * MROW * NCOL * MCOL C C -------------------------------------------- C ... FORM A(IBLK,JBLK) = A(IBLK,JBLK) C - W(IBLK) * Y(JBLK)**T C -------------------------------------------- C C... CALL HSMMPG ( IOPT, MROW, NCOL, MCOL, C... 1 WORK(IWM1+IROW), 1, MAXROW, C... 2 WORK(IYM1+JCOL), MAXROW, 1, C... 3 WORK(IR), 1, BLKSIZ, IERR ) C JPTR1 = IYM1+JCOL JPTR2 = IR DO 240 J = 1, MCOL CALL SGEMV ( 'N', MROW, NCOL, -1.0, 1 WORK(IWM1+IROW), MAXROW, 2 WORK(JPTR1), MAXROW, 3 1.0, WORK(JPTR2), 1 ) JPTR1 = JPTR1 + 1 JPTR2 = JPTR2 + BLKSIZ 240 CONTINUE C OPS = OPS + 2 * MROW * NCOL * MCOL C C ----------------------------------- C ... WRITE OUT MODIFIED A(IBLK,JBLK) C ----------------------------------- C WRITE ( UNIT1, REC=IREC, ERR=8100 ) 1 ( WORK(IRM1+J), J = 1, BLKSIZ**2 ) C 300 CONTINUE C C --------------------------------------- C ... END OF WORK ON CURRENT BLOCK COLUMN C --------------------------------------- C 400 CONTINUE C C----------------------------------------------------------------------- C C ------------------------------------------------------ C ... READ IN BANDED FORM, COMPUTE IT'S NORM AND COMPUTE C ITS EIGENVALUES AND EIGENVECTORS C ------------------------------------------------------ C MBAND = BLKSIZ + 1 LDA = NEQNS IF ( MOD ( LDA, 2 ) .EQ. 0 ) LDA = LDA + 1 C IS = 1 IR = IS + LDA * MBAND C CALL HSSXE4 ( NUMBLK, BLKSIZ, LDA, NEQNS, MBAND, WORK(IS), 1 WORK(IR), BLKSIZ**2, UNIT1, BNORM, IERR ) IF ( IERR .NE. 0 ) GO TO 8100 C C ----------------------------------------------------- C ... REDUCE TO TRIDIAGONAL FORM USING GIVENS ROTATIONS C ----------------------------------------------------- C CALL BANDR1 ( LDA, NEQNS, MBAND, WORK(IS), EIGVAL, 1 WORK(IR), WORK(IR+NEQNS), .FALSE., DUMMY ) C C ----------------------------------------------- C ... COMPUTE EIGENVALUES OF THE TRIDIAGONAL FORM C ----------------------------------------------- C CALL TQLRAT ( NEQNS, EIGVAL, WORK(IR+NEQNS), IERR ) C GO TO 9000 C C----------------------------------------------------------------------- C C --------------- C ... ERROR TRAPS C --------------- C C ---------------------- C ... I/O ERROR ON UNIT1 C ---------------------- C 8100 IERR = -5 GO TO 9000 C C ---------------------- C ... I/O ERROR ON UNIT2 C ---------------------- C 8200 IERR = -6 GO TO 9000 C C----------------------------------------------------------------------- C C ----------------- C ... END OF HSSXE2 C ----------------- C 9000 RETURN END SUBROUTINE HSSXE3 ( IREC, LDA, BLKSIZ, NROW, NCOL, 1 A, TEMP, NTEMP, LUNIT, 2 IERR ) C C---------------------------------------------------------------------- C C HSSXE3 IS A LOW LEVEL SUBROUTINE FOR HSSXE2 WHICH READS RECORD C NO. IREC FROM C THE LOGICAL I/O UNIT LUNIT. THE MATRIX BLOCK IS CONTAINED IN THE C ARRAY A WHOSE LEADING DIMENSION IS LDA AND WHOSE SIZE IS NROW C BY NCOL. IERR IS AN ERROR RETURN FLAG. C C WRITTEN BY ROGER G GRIMES, BOEING COMPUTER SERVICE, BELLEVUE, C WASHINGTON, FEBRUARY, 1987. C C---------------------------------------------------------------------- C C ------------------------ C ... VARIABLE DECLARATION C ------------------------ C IMPLICIT COMPLEX ( A - Z ) C INTEGER IREC, LDA, BLKSIZ, NROW, NCOL, 1 NTEMP, LUNIT, IERR C REAL A(LDA,NCOL), TEMP(NTEMP) C C ------------------------------ C ... LOCAL VARIABLE DECLARATION C ------------------------------ C INTEGER ICOL, ISTART C C -------------------- C ... SUBPROGRAMS USED C -------------------- C EXTERNAL SCOPY C C---------------------------------------------------------------------- C C ---------------- C ... READ IN TEMP C ---------------- C READ ( LUNIT, REC=IREC, ERR=8000 ) TEMP C C -------------------- C ... COPY TEMP INTO A C -------------------- C ISTART = 1 C DO 100 ICOL = 1, NCOL CALL SCOPY ( NROW, TEMP(ISTART), 1, A(1,ICOL), 1 ) ISTART = ISTART + BLKSIZ 100 CONTINUE C GO TO 9000 C C----------------------------------------------------------------------- C C ------------------ C ... I/O ERROR TRAP C ------------------ C 8000 IERR = -1 C C----------------------------------------------------------------------- C C ---------------- C ... END OF HSSXE3 C ---------------- C 9000 RETURN END SUBROUTINE HSSXE4 ( NUMBLK, BLKSIZ, LDA, NEQNS, MBAND, 1 ABAND, TEMP, NTEMP, 2 LUNIT, BNORM, IERR ) C C---------------------------------------------------------------------- C C HSSXE4 IS A LOW LEVEL SUBROUTINE FOR HSSXE2 WHICH READS THE C BANDED FORM OF THE REDUCED MATRIX AND STORES C IT INTO THE ARRAY ABAND. THE REDUCED MATRIX AND THE UNITARY C TRANSFORMATIONS HAVE OVERWRITTEN THE ORIGINAL MATRIX ON C THE LOGICAL I/O UNIT LUNIT. C C WRITTEN BY ROGER G GRIMES, BOEING COMPUTER SERVICE, BELLEVUE, C WASHINGTON, FEBRUARY, 1987. C C---------------------------------------------------------------------- C C ------------------------ C ... VARIABLE DECLARATION C ------------------------ C IMPLICIT COMPLEX ( A - Z ) C INTEGER NUMBLK, BLKSIZ, LDA, NEQNS, MBAND, 1 NTEMP, LUNIT, IERR C REAL ABAND(LDA,MBAND), TEMP(NTEMP), BNORM C C ------------------------------ C ... LOCAL VARIABLE DECLARATION C ------------------------------ C INTEGER IBLK, IDIAG, IREC, ISTART, JBLK, 1 KBLK, KCOL, KROW, L, LL, 2 MROW, RECORD C C -------------------- C ... SUBPROGRAMS USED C -------------------- C REAL SASUM C EXTERNAL SASUM, SCOPY C C ---------------------- C ... STATEMENT FUNCTION C ---------------------- C RECORD(IBLK,JBLK) = ( (JBLK - 1) * (2*NUMBLK + 2 - JBLK) ) / 2 1 + ( IBLK - JBLK + 1 ) C C---------------------------------------------------------------------- C C ------------------------- C ... FOR EACH BLOCK COLUMN C ------------------------- C DO 400 KBLK = 1, NUMBLK C C ------------------------------ C ... READ IN THE DIAGONAL BLOCK C ------------------------------ C IREC = RECORD ( KBLK, KBLK ) READ ( LUNIT, REC=IREC, ERR=8000 ) TEMP C C ---------------------------------- C ... MOVE LOWER TRIANGLE INTO ABAND C ---------------------------------- C KROW = ( KBLK-1 ) * BLKSIZ MROW = MIN ( BLKSIZ, NEQNS - KROW ) KCOL = BLKSIZ + 2 L = MROW + 1 C DO 200 IDIAG = 1, MROW L = L - 1 KROW = KROW + 1 KCOL = KCOL - 1 CALL SCOPY ( L, TEMP(IDIAG), BLKSIZ+1, 1 ABAND(KROW,KCOL), 1 ) 200 CONTINUE C IF ( KBLK .EQ. NUMBLK ) GO TO 400 C C ---------------------------------- C ... READ IN THE SUB-DIAGONAL BLOCK C ---------------------------------- C IREC = IREC + 1 READ ( LUNIT, REC=IREC, ERR=8000 ) TEMP C C ---------------------------------- C ... MOVE UPPER TRIANGLE INTO ABAND C ---------------------------------- C MROW = MIN ( BLKSIZ, NEQNS - KROW ) KROW = KROW + 1 LL = BLKSIZ + 1 ISTART = 1 - BLKSIZ DO 300 IDIAG = 1, BLKSIZ LL = LL - 1 L = MIN ( MROW, LL ) ISTART = ISTART + BLKSIZ CALL SCOPY ( L, TEMP(ISTART), BLKSIZ+1, 1 ABAND(KROW,IDIAG), 1 ) 300 CONTINUE C 400 CONTINUE C C --------------------------------------- C ... COMPUTE THE NORM OF THE BAND MATRIX C --------------------------------------- C BNORM = 0.0E0 C DO 500 L = 1, MBAND C LL = MBAND + 1 - L C BNORM = MAX ( BNORM, SASUM ( NEQNS-LL+1, ABAND(LL,L), 1 ) ) C 500 CONTINUE C GO TO 9000 C C----------------------------------------------------------------------- C C ------------------ C ... I/O ERROR TRAP C ------------------ C 8000 IERR = -1 C C----------------------------------------------------------------------- C C ----------------- C ... END OF HSSXE4 C ----------------- C 9000 RETURN END SUBROUTINE HSSXE5 ( IREC, LDA, BLKSIZ, NROW, NCOL, 1 A, TEMP, NTEMP, LUNIT, 2 IERR ) C C---------------------------------------------------------------------- C C HSSXE5 IS A LOW LEVEL SUBROUTINE FOR HSSXE2 WHICH WRITES RECORD C NO. IREC TO C THE LOGICAL I/O UNIT LUNIT. THE MATRIX BLOCK IS CONTAINED IN THE C ARRAY A WHOSE LEADING DIMENSION IS LDA AND WHOSE SIZE IS NROW C BY NCOL. IERR IS AN ERROR RETURN FLAG. C C WRITTEN BY ROGER G GRIMES, BOEING COMPUTER SERVICE, BELLEVUE, C WASHINGTON, JANUARY, 1987. C C---------------------------------------------------------------------- C C ------------------------ C ... VARIABLE DECLARATION C ------------------------ C IMPLICIT COMPLEX ( A - Z ) C INTEGER IREC, LDA, BLKSIZ, NROW, NCOL, 1 NTEMP, LUNIT, IERR C REAL A(LDA,NCOL), TEMP(NTEMP) C C ------------------------------ C ... LOCAL VARIABLE DECLARATION C ------------------------------ C INTEGER ICOL, ISTART C C -------------------- C ... SUBPROGRAMS USED C -------------------- C EXTERNAL SCOPY C C---------------------------------------------------------------------- C C -------------------- C ... COPY A INTO TEMP C -------------------- C ISTART = 1 C DO 100 ICOL = 1, NCOL CALL SCOPY ( NROW, A(1,ICOL), 1, TEMP(ISTART), 1 ) ISTART = ISTART + BLKSIZ 100 CONTINUE C C ------------------ C ... WRITE OUT TEMP C ------------------ C WRITE ( LUNIT, REC=IREC, ERR=8000 ) TEMP C GO TO 9000 C C----------------------------------------------------------------------- C C ------------------ C ... I/O ERROR TRAP C ------------------ C 8000 IERR = -1 C C----------------------------------------------------------------------- C C ----------------- C ... END OF HSSXE5 C ----------------- C 9000 RETURN END SUBROUTINE HSSXE6 ( X, LDX, N, P, R, Y, V ) C C----------------------------------------------------------------------- C C HSSXE6 IS A LOW LEVEL SUBROUTINE OF HSSXE2 WHICH COMPUTES THE QR C DECOMPOSITION OF X AND ACCUMULATES THE HOUSEHOLDER TRANSFORMATIONS C INTO THE W-Y REPRESENTATION. X IS OVERWRITTEN WITH W, Y AND R ARE C STORED IN THE ARRAYS Y AND R. C C WRITTEN BY ROGER G. GRIMES, BOEING COMPUTER SERVICES, BELLEVUE, C WASHINGTON, FEBRUARY, 1987. C C----------------------------------------------------------------------- C C ------------------------ C ... VARIABLE DECLARATION C ------------------------ C IMPLICIT NONE C INTEGER LDX, N, P C REAL X(LDX,P), R(P,P), Y(LDX,P), 1 V(P) C C ------------------- C ... LOCAL VARIABLES C ------------------- C INTEGER J, L, LEND C LOGICAL QNULL C REAL ALPHA, DIAG, NRMXL C REAL OPS COMMON / RGGOPS / OPS C C -------------------- C ... SUBPROGRAMS USED C -------------------- C REAL SDOT, SNRM2 C EXTERNAL SAXPY, SCOPY, SDOT, SNRM2, SSCAL, 1 SRSCL C C----------------------------------------------------------------------- C C ------------------------ C ... FOR EACH COLUMN OF X C ------------------------ C LEND = MIN ( N, P ) C DO 600 L = 1, LEND C NRMXL = -X(L,L) QNULL = .FALSE. C IF ( L .EQ. N ) THEN QNULL = .TRUE. GO TO 500 END IF C C ---------------------------------------------------------- C ... COMPUTE THE HOUSEHOLDER TRANSFORMATION TO ZERO OUT THE C L-TH COLUMN OF X. C ---------------------------------------------------------- C DIAG = X(L,L) CALL SSCAL ( N, 0.0E0, Y(1,L), 1 ) C NRMXL = SNRM2 ( N-L+1, X(L,L), 1 ) C OPS = OPS + 2 * ( N - L + 1 ) C IF ( NRMXL .EQ. 0.0E0 ) THEN QNULL = .TRUE. GO TO 500 END IF C IF ( DIAG .NE. 0.0E0 ) NRMXL = SIGN ( NRMXL, DIAG ) C CALL SRSCL ( N-L+1, NRMXL, X(L,L), 1 ) C X(L,L) = 1.0E0 + X(L,L) C OPS = OPS + N - L + 2 C C ---------------------------------------------------- C ... APPLY THE TRANSFORMATION TO COLUMNS L+1 TO P AND C ROWS P TO N OF X (I.E. THE REMAINDER OF X). C ---------------------------------------------------- C ALPHA = X(L,L) C DO 100 J = L+1, P V(J) = SDOT(N-L+1, X(L,L), 1, X(L,J), 1 ) 100 CONTINUE C CALL SRSCL ( P-L, -ALPHA, V(L+1), 1 ) C DO 200 J = L+1, P CALL SAXPY ( N-L+1, V(J), X(L,L), 1, X(L,J), 1 ) 200 CONTINUE C DO 300 J = 1, L-1 V(J) = SDOT ( N-L+1, X(L,L), 1, Y(L,J), 1 ) 300 CONTINUE C CALL SRSCL ( L-1, -ALPHA, V, 1 ) C DO 400 J = 1, L-1 CALL SAXPY ( N-L+1, V(J), X(L,L), 1, Y(L,J), 1 ) 400 CONTINUE C OPS = OPS + 4 * ( P-1 ) * ( N-L+1 ) + P-1 C 500 CALL SCOPY ( L-1, X(1,L), 1, R(1,L), 1 ) R(L,L) = -NRMXL CALL SSCAL ( P-L, 0.0E0, R(L+1,L), 1 ) C CALL SSCAL ( L-1, 0.0E0, X(1,L), 1 ) C CALL SCOPY ( N-L+1, X(L,L), 1, Y(L,L), 1 ) C IF ( QNULL ) THEN CALL SSCAL ( N-L+1, 0.0, X(L,L), 1 ) ELSE CALL SRSCL ( N-L+1, ALPHA, X(L,L), 1 ) END IF C OPS = OPS + N - L + 1 C 600 CONTINUE C C ---------------------------------------------- C ... IF ( N .LT. P ) COPY REMAINDER OF X INTO R C ---------------------------------------------- C IF ( N .LT. P ) THEN C DO 700 L = N+1, P CALL SCOPY ( N, X(1,L), 1, R(1,L), 1 ) 700 CONTINUE C END IF C C----------------------------------------------------------------------- C C ----------------- C ... END OF HSSXE6 C ----------------- C RETURN END SUBROUTINE HSSXE7 ( I, NEQNS, MB, UNIT1, UNIT2, 1 UNIT3, NVEC, BNORM, EVAL, BANDMX, 2 EVEC, WORK, IERR ) C C----------------------------------------------------------------------- C C HSSXE7 IS A LOW LEVEL SUBROUTINE FOR HSSXE1. IT COMPUTES A BLOCK C OF EIGENVECTORS FOR A SPECIFIED SET OF EIGENVALUES. C C IF IT IS THE FIRST SET TO BE COMPUTED (I=1) THEN HSSXE7 READS THE C BANDED (REDUCED) FORM OF A FROM UNIT1 AND STORES IT INTO ARRAY C BANDMX. IT THEN INVOKES EISPACK SUBROUTINE BANDV TO COMPUTE C THE SPECIFIED EIGENVECTORS OF THE BANDED MATRIX. FINALLY THE C W-Y REPRESENTATION OF THE BLOCK HOUSEHOLDER TRANSFORMATIONS C USED TO REDUCE A TO BANDED FORM ARE BROUGHT IN AND THE C EIGENVECTORS ARE BACK TRANSFORMED TO THOSE OF THE DENSE MATRIX. C C HSSXE7 ASSUMES THAT THE BAND MATRIX HAS BEEN WRITTEN TO I/O UNIT C UNIT1 IN BLOCK FORM BY HSSXEV. FURTHERMORE IT ASSUMES THAT THE C W-Y REPRESENTATION IS STORED ON UNIT2. C C WRITTEN BY ROGER G. GRIMES, BOEING COMPUTER SERVICES, BELLEVUE, C WASHINGTON, MARCH, 1987. C C----------------------------------------------------------------------- C C ------------------------ C ... VARIABLE DECLARATION C ------------------------ C IMPLICIT COMPLEX ( A - Z ) C INTEGER I, NEQNS, MB, UNIT1, UNIT2, 1 UNIT3, NVEC, IERR C REAL BNORM C REAL EVAL(NVEC), BANDMX(NEQNS,MB), 1 EVEC(NEQNS,NVEC), WORK(*) C C ------------------- C ... LOCAL VARIABLES C ------------------- C INTEGER BLKSIZ, IBLK, IREC, IREC2, IROW, 1 IT1, IT2, J, JBLK, JPTR, 2 JROW, K, KBLK, MAXROW, MROW, 3 NROW, NUMBLK, RECORD C LOGICAL QOPEN C REAL V C REAL OPS COMMON / RGGOPS / OPS C C -------------------- C ... SUBPROGRAMS USED C -------------------- C REAL SNRM2 C EXTERNAL BANDVM, HSSXE3, HSSXE4, SNRM2, SGEMV, 1 SSCAL C C ---------------------- C ... STATEMENT FUNCTION C ---------------------- C RECORD(IBLK,JBLK) = ( (JBLK - 1) * (2*NUMBLK + 2 - JBLK) ) / 2 1 + ( IBLK - JBLK + 1 ) C C----------------------------------------------------------------------- C BLKSIZ = MB - 1 NUMBLK = ( NEQNS + BLKSIZ - 1 ) / BLKSIZ MAXROW = NEQNS - BLKSIZ C C -------------------------------------------------------- C ... IF FIRST CLUSTER PROCESSED THEN READ IN BAND MATRIX. C -------------------------------------------------------- C IF ( I .EQ. 1 ) THEN C C ----------------------- C ... READ IN BANDED FORM C ----------------------- C CALL HSSXE4 ( NUMBLK, BLKSIZ, NEQNS, NEQNS, MB, BANDMX, 1 WORK, BLKSIZ**2, UNIT1, BNORM, IERR ) IF ( IERR .NE. 0 ) GO TO 8100 C C ----------------------------------- C ... INQUIRE ON THE STATUS OF UNIT3. C ----------------------------------- C INQUIRE ( UNIT=UNIT3, ERR=8300, OPENED=QOPEN ) C IF ( QOPEN ) THEN C C --------------------------------------- C ... UNIT3 IS ALREADY OPENED. CLOSE IT. C --------------------------------------- C CLOSE ( UNIT=UNIT3, ERR=8300 ) C END IF C C --------------- C ... OPEN UNIT3. C --------------- C OPEN ( UNIT=UNIT3, ERR=8300, ACCESS='DIRECT', 1 STATUS='SCRATCH', RECL=8*NEQNS ) C END IF C C----------------------------------------------------------------------- C C ----------------------------------------------------------- C ... COMPUTE THE CLUSTER OF EIGENVECTORS FOR THE BAND MATRIX C ----------------------------------------------------------- C IT1 = NEQNS * MB + 1 C CALL BANDVM ( NEQNS, NEQNS, MB, BANDMX, BNORM, NVEC, EVAL, 1 EVEC, IERR, WORK, WORK(IT1) ) C C -------------------------------------------------- C ... BACK TRANSFORM THE EIGENVECTORS OF THE BANDED C MATRIX TO THE EIGENVECTORS OF THE DENSE MATRIX C -------------------------------------------------- C IT1 = MAXROW * BLKSIZ + 1 IT2 = IT1 + BLKSIZ **2 C C ------------------ C ... FOR EACH BLOCK C ------------------ C DO 800 KBLK = NUMBLK-1, 1, -1 C C ------------- C ... READ IN Y C ------------- C IREC2 = 2 * ( RECORD ( KBLK, KBLK ) - KBLK ) IROW = 1 - BLKSIZ C JROW = KBLK * BLKSIZ + 1 NROW = NEQNS - JROW + 1 C DO 700 IBLK = KBLK + 1, NUMBLK C IROW = IROW + BLKSIZ MROW = MIN ( NROW-IROW+1, BLKSIZ ) C IREC2 = IREC2 + 1 CALL HSSXE3 ( IREC2, MAXROW, BLKSIZ, MROW, BLKSIZ, 1 WORK(IROW), WORK(IT2), BLKSIZ**2, 2 UNIT2, IERR ) IF ( IERR .NE. 0 ) GO TO 8200 C 700 CONTINUE C C --------------------- C ... FORM T1 = Y' * X C --------------------- C C... CALL HSMMPG ( 11, BLKSIZ, NROW, NVEC, WORK, 1, MAXROW, C... 1 EVEC(JROW,1), 1, NEQNS, WORK(IT1), 1, BLKSIZ, C... 2 IERR ) C JPTR = IT1 DO 710 J = 1, NVEC CALL SGEMV ( 'T', NROW, BLKSIZ, 1.0, WORK, MAXROW, 1 EVEC(JROW,J), 1, 0.0, WORK(JPTR), 1 ) JPTR = JPTR + BLKSIZ 710 CONTINUE C OPS = OPS + 2 * BLKSIZ * NROW * NVEC C C ------------- C ... READ IN W C ------------- C IROW = 1 - BLKSIZ C JROW = KBLK * BLKSIZ + 1 NROW = NEQNS - JROW + 1 C DO 750 IBLK = KBLK + 1, NUMBLK C IROW = IROW + BLKSIZ MROW = MIN ( NROW-IROW+1, BLKSIZ ) C IREC2 = IREC2 + 1 CALL HSSXE3 ( IREC2, MAXROW, BLKSIZ, MROW, BLKSIZ, 1 WORK(IROW), WORK(IT2), BLKSIZ**2, 2 UNIT2, IERR ) IF ( IERR .NE. 0 ) GO TO 8200 C 750 CONTINUE C C ----------------------- C ... FORM X = X - W * T1 C ----------------------- C C... CALL HSMMPG ( 3, NROW, BLKSIZ, NVEC, WORK, 1, MAXROW, C... 1 WORK(IT1), 1, BLKSIZ, EVEC(JROW,1), 1, NEQNS, C... 2 IERR ) C JPTR = IT1 DO 790 J = 1, NVEC CALL SGEMV ( 'N', NROW, BLKSIZ, -1.0, WORK, MAXROW, 1 WORK(JPTR), 1, 1.0, EVEC(JROW,J), 1 ) JPTR = JPTR + BLKSIZ 790 CONTINUE C OPS = OPS + 2 * NROW * BLKSIZ * NVEC C 800 CONTINUE C C -------------------------------------------------------- C ... NORMALIZE AND WRITE OUT THE TRANSFORMED EIGENVECTORS C -------------------------------------------------------- C IREC = I - 1 DO 900 J = 1, NVEC C V = 1.0E0 / SNRM2 ( NEQNS, EVEC(1,J), 1 ) CALL SSCAL ( NEQNS, V, EVEC(1,J), 1 ) C IREC = IREC + 1 WRITE ( UNIT3, REC=IREC, ERR=8300 ) 1 ( EVEC(K,J), K = 1, NEQNS ) C 900 CONTINUE C OPS = OPS + 3 * NVEC * NEQNS C GO TO 9000 C C----------------------------------------------------------------------- C C --------------- C ... ERROR TRAPS C --------------- C C ---------------------- C ... I/O ERROR ON UNIT1 C ---------------------- C 8100 IERR = -6 GO TO 9000 C C ---------------------- C ... I/O ERROR ON UNIT2 C ---------------------- C 8200 IERR = -7 GO TO 9000 C C ---------------------- C ... I/O ERROR ON UNIT3 C ---------------------- C 8300 IERR = -8 GO TO 9000 C C----------------------------------------------------------------------- C C ----------------- C ... END OF HSSXE7 C ----------------- C 9000 RETURN END SUBROUTINE BANDR1 (NM,N,MB,A,D,E,E2) C INTEGER J,K,N,R,I1,I2,J1,KR,MB,MR,M1,NM,MAXL,MAXR, 1 LL, ROWBGN REAL A(NM,MB),D(N),E(N),E2(N) REAL A11, A12, A21, A22, C, G, S C REAL OPS COMMON / RGGOPS / OPS C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BANDRD, C NUM. MATH. 12, 231-241(1968) BY SCHWARZ. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 273-283(1971). C C THIS SUBROUTINE REDUCES A REAL SYMMETRIC BAND MATRIX C TO A SYMMETRIC TRIDIAGONAL MATRIX USING ORTHOGONAL SIMILARITY C 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 MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE C NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL C DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE C LOWER TRIANGLE OF THE MATRIX. C C A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT C MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL C IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, C ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE C SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY C ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN. C CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. C C ON OUTPUT C C A HAS BEEN DESTROYED, EXCEPT FOR ITS LAST TWO COLUMNS WHICH C CONTAIN A COPY OF THE TRIDIAGONAL MATRIX. 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 E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. 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 IS MODIFIED TO USE SROTG AND SROT FROM BLAS C TO IMPROVE PERFORMANCE ON THE CRAY X-MP. MODIFICATIONS C BY ROGER G. GRIMES, BOEING COMPUTER SERVICES, BELLEVUE, C WASHINGTON, FEBRUARY, 1987. C C ------------------------------------------------------------------ C IF ( N .LE. 0 .OR. MB .LE. 0 ) RETURN C IF ( MB .EQ. 1 ) GO TO 900 IF ( MB .EQ. 2 ) GO TO 800 C M1 = MB - 1 C DO 700 K = 1, N-2 MAXR = MIN (M1,N-K) C .......... FOR R=MAXR STEP -1 UNTIL 2 DO -- .......... DO 600 R = MAXR, 2, -1 KR = K + R MR = MB - R G = A(KR,MR) ROWBGN = MR + 1 C DO 500 J = KR, N, M1 J1 = J - 1 IF (G .EQ. 0.0E0) GO TO 600 C C ... COMPUTE THE GIVENS ROTATION TO ANNIHALATE C G USING ELEMENT A(J1,ROWBGN). C CALL SROTG ( A(J1,ROWBGN), G, C, S ) C C ... APPLY THE ROTATION TO THE DIAGONAL BLOCK C A11 = C * A(J1,MB) + S * A(J,M1) A12 = C * A(J ,M1) + S * A(J,MB) A21 = -S * A(J1,MB) + C * A(J,M1) A22 = -S * A(J ,M1) + C * A(J,MB) C A(J1,MB) = C * A11 + S * A12 A(J ,M1) = C * A21 + S * A22 A(J ,MB) = -S * A21 + C * A22 C C ... APPLY THE ROTATION TO ROWS J1 AND J C LL = MB - ROWBGN - 1 CALL SROT ( LL, A(J1,ROWBGN+1), NM, 1 A(J, ROWBGN ), NM, C, S ) C OPS = OPS + 6 * LL + 31 C C ... APPLY THE ROTATION TO COLUMNS J1 AND J C IF (J .EQ. N) GO TO 500 C MAXL = MIN ( MB-2, N-J ) I1 = J + MAXL I2 = M1 - MAXL C CALL SROT ( MAXL, A(I1,I2 ), NM-1, 1 A(I1,I2+1), NM-1, C, S ) C OPS = OPS + 6 * MAXL C C ... COMPUTE THE INTRODUCED SPURIOUS ELEMENT C I1 = I1 + 1 IF ( I1 .LE. N ) THEN G = S * A(I1,1) A(I1,1) = C * A(I1,1) OPS = OPS + 2 END IF C ROWBGN = 1 C 500 CONTINUE C 600 CONTINUE C 700 CONTINUE C C ... SAVE THE TRIDIAGONAL FORM IN ARRAYS D AND E C 800 D(1) = A(1,MB) E (1) = 0.0E0 E2(1) = 0.0E0 C DO 830 J = 2, N D (J) = A(J,MB ) E (J) = A(J,MB-1) E2(J) = A(J,MB-1) ** 2 830 CONTINUE C OPS = OPS + ( N - 1 ) C RETURN C C ... HANDLE DIAGONAL MATRIX CASE C 900 DO 950 J = 1, N D(J) = A(J,MB) E(J) = 0.0E0 E2(J) = 0.0E0 950 CONTINUE C 1001 RETURN END SUBROUTINE BANDVM(NM,N,MB,A,NORM,M,W,Z,IERR,ASHIFT,RV6) C********************************************* C * * C * COPYRIGHT 1987 THE BOEING COMPANY * C * * C * ALL RIGHTS RESERVED * C * * C ******************************************** C INTEGER I,J,K,M,N,R,JJ,M1,NM,ITS,MB,IERR,GROUP C REAL A(NM,MB),W(M),Z(NM,M),ASHIFT(MB,N),RV6(N) REAL XU,X0,X1,EPS2,EPS3,EPS4,NORM,MACEPS REAL SDOT, SMACH, SNRM2 C REAL OPS COMMON / RGGOPS / OPS C C THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A REAL SYMMETRIC C BAND MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, USING INVERSE C ITERATION. 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 MB IS THE NUMBER OF COLUMNS OF THE ARRAY A USED TO STORE THE C BAND MATRIX, I.E. THE NUMBER OF SUBDIAGONALS PLUS 1 (FOR C THE PRINCIPAL DIAGONAL). C C A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT C MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL C IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, C ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE C SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY C ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF COLUMN MB. C CONTENTS OF A NOT PART OF THE MATRIX ARE ARBITRARY. C C M IS THE NUMBER OF SPECIFIED EIGENVALUES. C C W CONTAINS THE M EIGENVALUES IN ASCENDING ORDER. C C WORKING STORAGE C C ASHIFT IS USED TO HOLD A - LAMBDA*I FOR EACH APPLICATION C OF INVERSE ITERATION. C RV6 IS A VECTOR OF LENGTH N USED IN INVERSE ITERATION. C C ON OUTPUT C C A AND W ARE UNALTERED. C C Z CONTAINS THE ASSOCIATED SET OF ORTHOGONAL EIGENVECTORS. C ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO. C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH C EIGENVALUE FAILS TO CONVERGE. C 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 MODIFIED HEAVILY BY ROGER G. GRIMES, BOEING COMPUTER C SERVICES, MARCH, 1987. C C ------------------------------------------------------------------ C MACEPS = SMACH(1) IERR = 0 IF (M .EQ. 0) GO TO 1001 M1 = MB - 1 X0 = W(1) C C .......... EPS2 IS THE CRITERION FOR GROUPING, C EPS3 REPLACES ZERO PIVOTS AND EQUAL C ROOTS ARE MODIFIED BY EPS3, C EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .......... C IF (NORM .EQ. 0.0E0) NORM = 1.0E0 EPS2 = 1.0E-3 * NORM EPS3 = MACEPS * NORM EPS4 = SQRT ( REAL ( N ) ) * EPS3 C C .......... FIND VECTORS BY INVERSE ITERATION .......... C DO 920 R = 1, M ITS = 1 X1 = W(R) C C .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... C IF (ABS(X1-X0) .GE. EPS2 .OR. R.EQ.1) THEN GROUP = 0 ELSE GROUP = GROUP + 1 IF ( X1 .LE. X0 ) X1 = X0 + EPS3 END IF C C .......... TRANSPOSE A INTO ASHIFT AND SUBTRACT EIGENVALUE C FROM MAIN DIAGONAL. ALSO INITIALIZE VECTOR ........ C DO 150 I = 1, M1 CALL SCOPY ( N, A(1,I), 1, ASHIFT(I,1), MB ) 150 CONTINUE C DO 200 I = 1, N ASHIFT(MB,I) = A(I,MB) - X1 RV6(I) = EPS4 200 CONTINUE C OPS = OPS + N C IF (M1 .EQ. 0) GO TO 600 C C .......... FACTOR ASHIFT INTO L*D*L' .......... C DO 320 J = 2, N LM = MAX ( J-M1, 1 ) LN = MAX ( MB+1-J, 1 ) LEN = MB - LN CALL STRSV ( 'UPPER', 'TRANSPOSED', 'UNIT', 1 LEN, ASHIFT(MB,LM), M1, ASHIFT(LN,J), 1 ) C OPS = OPS + LEN * ( LEN-1 ) C K = J - LEN - 1 CDIR$ IVDEP DO 300 I = LN, M1 K = K + 1 ASHIFT(I,J) = ASHIFT(I,J) / ASHIFT(MB,K) 300 CONTINUE C S = ASHIFT(MB,J) K = J - LEN - 1 DO 310 I = LN, M1 K = K + 1 S = S - ASHIFT(MB,K) * ASHIFT(I,J)**2 310 CONTINUE C OPS = OPS + 4 * ( M1 - LN + 1 ) C IF ( ABS(S) .LT. EPS3 ) S = SIGN(EPS3,S) ASHIFT(MB,J) = S C 320 CONTINUE C C .......... SOLVE WITH FACTORS OF ASHIFT .......... C 600 DO 620 K = 1, N-1 S = RV6(K) DO 610 J = 1, MIN ( M1, N-K ) RV6(K+J) = RV6(K+J) - S * ASHIFT(MB-J,K+J) 610 CONTINUE OPS = OPS + 2 * MIN ( M1, N-K ) 620 CONTINUE C DO 630 K = 1, N RV6(K) = RV6(K) / ASHIFT(MB,K) 630 CONTINUE C OPS = OPS + N C DO 650 K = N, 2, -1 S = RV6(K) CDIR$ IVDEP DO 640 J = 1, MIN ( K-1, M1 ) RV6(K-J) = RV6(K-J) - S * ASHIFT(MB-J,K) 640 CONTINUE OPS = OPS + 2 * MIN ( K-1, M1 ) 650 CONTINUE C C .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS C MEMBERS OF GROUP .......... C DO 680 JJ = 1, GROUP J = R - GROUP - 1 + JJ C XU = SDOT ( N, RV6, 1, Z(1,J), 1 ) C DO 660 I = 1, N 660 RV6(I) = RV6(I) - XU * Z(I,J) C OPS = OPS + 4 * N C 680 CONTINUE C C .......... COMPUTE NORM OF VECTOR AND TEST FOR C CONVERGENCE ......... C 700 NORM = SASUM ( N, RV6, 1 ) OPS = OPS + N C IF (NORM .GE. 0.1E0) GO TO 840 C C .......... EIGENVECTOR DID NOT CONVERGE. SCALE VECTOR C AND PERFORM ANOTHER ITERATION .......... C IF (ITS .GE. N) GO TO 830 ITS = ITS + 1 C IF ( NORM .EQ. 0.0E0 ) THEN C RV6(ITS) = EPS4 C ELSE C XU = EPS4 / NORM C DO 760 I = 1, N RV6(I) = RV6(I) * XU 760 CONTINUE C OPS = OPS + N + 1 C END IF C GO TO 600 C C .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... C 830 IERR = -R XU = 0.0E0 GO TO 870 C C .......... NORMALIZE SO THAT SUM OF SQUARES IS C 1 AND EXPAND TO FULL ORDER .......... C 840 XU = 1.0E0 / SNRM2 ( N, RV6, 1 ) OPS = OPS + 2 * N C 870 DO 900 I = 1, N 900 Z(I,R) = RV6(I) * XU OPS = OPS + N C X0 = X1 920 CONTINUE C 1001 RETURN END SUBROUTINE TQLRAT(N,D,E2,IERR) C INTEGER I,J,L,M,N,II,L1,MML,IERR REAL D(N),E2(N) REAL B,C,F,G,H,P,R,S,T,MACEPS,SMACH C REAL OPS COMMON / RGGOPS / OPS C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT, C ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH. C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC C TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE C INPUT MATRIX IN ITS LAST N-1 POSITIONS. E2(1) IS ARBITRARY. C C ON OUTPUT C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE C THE SMALLEST EIGENVALUES. C C E2 HAS BEEN DESTROYED. 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 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 MACEPS = SMACH(1) IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E2(I-1) = E2(I) C F = 0.0E0 T = 0.0E0 E2(N) = 0.0E0 C DO 290 L = 1, N J = 0 H = ABS(D(L)) + SQRT(E2(L)) OPS = OPS + 2 IF (T .GT. H) GO TO 105 T = H B = MACEPS * ABS(T) C = B * B OPS = OPS + 3 C .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT .......... 105 DO 110 M = L, N IF (E2(M) .LE. C) GO TO 120 C .......... E2(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 210 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C .......... FORM SHIFT .......... L1 = L + 1 S = SQRT(E2(L)) G = D(L) P = (D(L1) - G) / (2.0E0 * S) R = SQRT(P**2+1.0E0) D(L) = S / (P + SIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C OPS = OPS + 11 + N - L C C .......... RATIONAL QL TRANSFORMATION .......... G = D(M) IF (G .EQ. 0.0E0) G = B H = G S = 0.0E0 MML = M - L C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... DO 200 II = 1, MML I = M - II P = G * H R = P + E2(I) E2(I+1) = S * R S = E2(I) / R D(I+1) = H + S * (H + D(I)) G = D(I) - E2(I) / G IF (G .EQ. 0.0E0) G = B H = G * P / R 200 CONTINUE C E2(L) = S * G D(L) = H C OPS = OPS + 10 * MML + 1 C .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST .......... IF (H .EQ. 0.0E0) GO TO 210 IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210 E2(L) = H * E2(L) IF (E2(L) .NE. 0.0E0) GO TO 130 210 P = D(L) + F C .......... ORDER EIGENVALUES .......... IF (L .EQ. 1) GO TO 250 C .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... DO 230 II = 2, L I = L + 2 - II IF (P .GE. D(I-1)) GO TO 270 D(I) = D(I-1) 230 CONTINUE C 250 I = 1 270 D(I) = P 290 CONTINUE C GO TO 1001 C .......... SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS .......... 1000 IERR = L 1001 RETURN END REAL FUNCTION SASUM(N,SX,INCX) C C TAKES THE SUM OF THE ABSOLUTE VALUES. C REAL SX(1),STEMP INTEGER I,INCX,N,NINCX C SASUM = 0.0E0 STEMP = 0.0E0 IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX STEMP = STEMP + ABS(SX(I)) 10 CONTINUE SASUM = STEMP RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DO 30 I = 1,N STEMP = STEMP + ABS(SX(I)) 30 CONTINUE SASUM = STEMP RETURN END SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY) C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C REAL SX(1),SY(1),SA INTEGER I,INCX,INCY,IX,IY,N C IF(N.LE.0)RETURN IF (SA .EQ. 0.0) RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N SY(IY) = SY(IY) + SA*SX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C 20 DO 30 I = 1,N SY(I) = SY(I) + SA*SX(I) 30 CONTINUE RETURN END SUBROUTINE SCOPY(N,SX,INCX,SY,INCY) C C COPIES A VECTOR, X, TO A VECTOR, Y. C REAL SX(1),SY(1) INTEGER I,INCX,INCY,IX,IY,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N SY(IY) = SX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C 20 DO 30 I = 1,N SY(I) = SX(I) 30 CONTINUE RETURN END REAL FUNCTION SDOT(N,SX,INCX,SY,INCY) C C FORMS THE DOT PRODUCT OF TWO VECTORS. C REAL SX(1),SY(1),STEMP INTEGER I,INCX,INCY,IX,IY,N C STEMP = 0.0E0 SDOT = 0.0E0 IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N STEMP = STEMP + SX(IX)*SY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE SDOT = STEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C 20 DO 30 I = 1,N STEMP = STEMP + SX(I)*SY(I) 30 CONTINUE SDOT = STEMP RETURN END REAL FUNCTION SNRM2 ( N, SX, INCX) INTEGER NEXT REAL SX(1), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE DATA ZERO, ONE /0.0E0, 1.0E0/ C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN SX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C C.L.LAWSON, 1978 JAN 08 C C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF SQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF SQRT(V) OVER ALL KNOWN MACHINES. C WHERE C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) C V = LARGEST NO. (OVERFLOW LIMIT) C C BRIEF OUTLINE OF ALGORITHM.. C C PHASE 1 SCANS ZERO COMPONENTS. C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. C C VALUES FOR CUTLO AND CUTHI.. C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE C UNIVAC AND DEC AT 2**(-103) C THUS CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. C THUS CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. C THUS CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / C IF(N .GT. 0) GO TO 10 SNRM2 = ZERO GO TO 300 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C BEGIN MAIN LOOP I = 1 20 GO TO NEXT,(30, 50, 70, 110) 30 IF( ABS(SX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF( SX(I) .EQ. ZERO) GO TO 200 IF( ABS(SX(I)) .GT. CUTLO) GO TO 85 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 105 C C PREPARE FOR PHASE 4. C 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / SX(I)) / SX(I) 105 XMAX = ABS(SX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF( ABS(SX(I)) .GT. CUTLO ) GO TO 75 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 110 IF( ABS(SX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / SX(I))**2 XMAX = ABS(SX(I)) GO TO 200 C 115 SUM = SUM + (SX(I)/XMAX)**2 GO TO 200 C C C PREPARE FOR PHASE 3. C 75 SUM = (SUM * XMAX) * XMAX C C C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C 85 HITEST = CUTHI/FLOAT( N ) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 95 J =I,NN,INCX IF(ABS(SX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + SX(J)**2 SNRM2 = SQRT( SUM ) GO TO 300 C 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C SNRM2 = XMAX * SQRT(SUM) 300 CONTINUE RETURN END SUBROUTINE SROT (N,SX,INCX,SY,INCY,C,S) C C APPLIES A PLANE ROTATION. C REAL SX(1),SY(1),STEMP,C,S INTEGER I,INCX,INCY,IX,IY,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL C TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N STEMP = C*SX(IX) + S*SY(IY) SY(IY) = C*SY(IY) - S*SX(IX) SX(IX) = STEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C 20 DO 30 I = 1,N STEMP = C*SX(I) + S*SY(I) SY(I) = C*SY(I) - S*SX(I) SX(I) = STEMP 30 CONTINUE RETURN END SUBROUTINE SROTG(SA,SB,C,S) C C CONSTRUCT GIVENS PLANE ROTATION. C REAL SA,SB,C,S,ROE,SCALE,R,Z C ROE = SB IF( ABS(SA) .GT. ABS(SB) ) ROE = SA SCALE = ABS(SA) + ABS(SB) IF( SCALE .NE. 0.0 ) GO TO 10 C = 1.0 S = 0.0 R = 0.0 GO TO 20 10 R = SCALE*SQRT((SA/SCALE)**2 + (SB/SCALE)**2) R = SIGN(1.0,ROE)*R C = SA/R S = SB/R 20 Z = 1.0 IF( ABS(SA) .GT. ABS(SB) ) Z = S IF( ABS(SB) .GE. ABS(SA) .AND. C .NE. 0.0 ) Z = 1.0/C SA = R SB = Z RETURN END SUBROUTINE SRSCL(N,SA,SX,INCX) C C SCALES A VECTOR BY THE RECIPROCAL OF A CONSTANT. C REAL SA,SX(1) INTEGER I,INCX,N,NINCX C IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX SX(I) = SX(I) / SA 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DO 30 I = 1,N SX(I) = SX(I) / SA 30 CONTINUE RETURN END SUBROUTINE SSCAL(N,SA,SX,INCX) C C SCALES A VECTOR BY A CONSTANT. C REAL SA,SX(1) INTEGER I,INCX,N,NINCX C IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX SX(I) = SA*SX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DO 30 I = 1,N SX(I) = SA*SX(I) 30 CONTINUE RETURN END REAL FUNCTION SMACH(JOB) INTEGER JOB C C SMACH COMPUTES MACHINE PARAMETERS OF FLOATING POINT C ARITHMETIC FOR USE IN TESTING ONLY. NOT REQUIRED BY C LINPACK PROPER. C C IF TROUBLE WITH AUTOMATIC COMPUTATION OF THESE QUANTITIES, C THEY CAN BE SET BY DIRECT ASSIGNMENT STATEMENTS. C ASSUME THE COMPUTER HAS C C B = BASE OF ARITHMETIC C T = NUMBER OF BASE B DIGITS C L = SMALLEST POSSIBLE EXPONENT C U = LARGEST POSSIBLE EXPONENT C C THEN C C EPS = B**(1-T) C TINY = 100.0*B**(-L+T) C HUGE = 0.01*B**(U-T) C C DMACH SAME AS SMACH EXCEPT T, L, U APPLY TO C DOUBLE PRECISION. C C CMACH SAME AS SMACH EXCEPT IF COMPLEX DIVISION C IS DONE BY C C 1/(X+I*Y) = (X-I*Y)/(X**2+Y**2) C C THEN C C TINY = SQRT(TINY) C HUGE = SQRT(HUGE) C C C JOB IS 1, 2 OR 3 FOR EPSILON, TINY AND HUGE, RESPECTIVELY. C C REAL EPS,TINY,HUGE,S C EPS = 1.0 10 EPS = EPS/2.0 S = 1.0 + EPS IF (S .GT. 1.0) GO TO 10 EPS = 2.0*EPS C S = 1.0 20 TINY = S S = S/16.0 IF (S*100. .NE. 0.0) GO TO 20 TINY = (TINY/EPS)*100.0 HUGE = 1.0/TINY C IF (JOB .EQ. 1) SMACH = EPS IF (JOB .EQ. 2) SMACH = TINY IF (JOB .EQ. 3) SMACH = HUGE RETURN END SUBROUTINE SGTHR ( NZ, Y, X, INDX ) C C ================================================================== C ================================================================== C ==== SGTHR -- REAL GATHER ==== C ================================================================== C ================================================================== C C PURPOSE C ------- C C SGTHR GATHERS THE SPECIFIED ELEMENTS FROM C A REAL VECTOR Y IN FULL STORAGE FORM C INTO C A REAL VECTOR X IN COMPRESSED FORM (X,INDX). C C ONLY THE ELEMENTS OF Y WHOSE INDICES ARE LISTED IN INDX C ARE REFERENCED. C C ARGUMENTS C --------- C C INPUT ... C C NZ INTEGER NUMBER OF ELEMENTS TO BE GATHERED INTO C COMPRESSED FORM. C Y REAL ARRAY, ON INPUT, WHICH CONTAINS THE C VECTOR Y IN FULL STORAGE FORM. ONLY C THE ELEMENTS CORRESPONDING TO THE INDICES C IN INDX WILL BE ACCESSED. C INDX INTEGER ARRAY CONTAINING THE INDICES OF THE C VALUES TO BE GATHERED INTO COMPRESSED FORM. C C OUTPUT ... C C X REAL ARRAY CONTAINING THE VALUES GATHERED INTO C THE COMPRESSED FORM. C C SPARSE BASIC LINEAR ALGEBRA SUBPROGRAM C C FORTRAN VERSION WRITTEN OCTOBER 1984 C ROGER G GRIMES, BOEING COMPUTER SERVICES C C ================================================================== C C ------------- C ... ARGUMENTS C ------------- C C INTEGER NZ, INDX (*) C REAL Y (*), X (*) C C ------------------- C ... LOCAL VARIABLES C ------------------- C INTEGER I C C ================================================================== C IF ( NZ .LE. 0 ) RETURN C DO 10 I = 1, NZ X(I) = Y(INDX(I)) 10 CONTINUE C RETURN END SUBROUTINE SGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. Scalar Arguments .. REAL ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. Array ARguments .. REAL A( LDA, * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * SGEMV 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 . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL 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 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 . * 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 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 20-July-1986. * Sven Hammarling, NAG Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. REAL ONE , ZERO PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) * .. Local Scalars .. REAL 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( 'SGEMV ', 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. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N 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 and set up the start points in X and Y * if the increments are not both unity. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN IF( BETA.NE.ONE )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 END IF ELSE 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 IF( BETA.NE.ONE )THEN 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. * IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN DO 60, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = ALPHA*X( J ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF 60 CONTINUE ELSE JX = KX 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. * IF( ( INCX.EQ.1 ).AND.( INCY.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( J ) = Y( J ) + ALPHA*TEMP 100 CONTINUE ELSE JY = KY 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 SGEMV . * END SUBROUTINE STRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. Scalar Arguments .. INTEGER INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. Array Arguments .. REAL A( LDA, * ), X( * ) * .. * * PURPOSE * ======= * * STRSV Solves one of the systems of equations * * A*X = B, or A'*X = B, * * where B and X are N element vectors and A is an N by N unit, or * non-unit, upper or lower triangular matrix. * * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. * * PARAMETERS * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' or 'u' A is an upper triangular matrix. * * UPLO = 'L' or 'l' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the equations to be solved as * follows: * * TRANS = 'N' or 'n' A*X = B. * * TRANS = 'T' or 't' A'*X = B. * * TRANS = 'C' or 'c' A'*X = B. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' or 'u' A is assumed to be unit triangular. * * DIAG = 'N' or 'n' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, N ). * Before entry with UPLO = 'U' or 'u', the leading N by N * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L' or 'l', the leading N by N * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U' or 'u', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least * MAX( 1, N ). * Unchanged on exit. * * X - REAL array of DIMENSION at least * ( 1 + ( N - 1 )*ABS( INCX ) ). * Before entry, the incremented array X must contain the N * element right-hand side vector B. On exit, X is * overwritten with the solution vector X. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. INCX must not be zero. * Unchanged on exit. * * * LEVEL 2 BLAS ROUTINE. * * -- Written on 20-July-1986. * Sven Hammarling, NAG Central Office. * Richard Hanson, Sandia National Labs. * * * .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) * .. Local Scalars .. INTEGER I, INFO, IX, J, JX, KX LOGICAL NOUNIT * .. 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( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = 1 ELSE IF ( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) ) THEN INFO = 2 ELSE IF ( .NOT.LSAME( DIAG, 'U' ).AND. $ .NOT.LSAME( DIAG, 'N' ) ) THEN INFO = 3 ELSE IF ( N.LT.0 ) THEN INFO = 4 ELSE IF ( LDA.LT.MAX(1,N) ) THEN INFO = 6 ELSE IF ( INCX.EQ.0 ) THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'STRSV ', INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * NOUNIT = LSAME( DIAG, 'N' ) * * Set up the start point in X if the increment is not unity. * This will be ( N - 1 )*INCX too small for descending loops. * IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through a. * IF( LSAME( TRANS, 'N' ) )THEN * * Form X := INV( A )*X. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 20, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) DO 10, I = J - 1, 1, -1 X( I ) = X( I ) - X( J )*A( I, J ) 10 CONTINUE END IF 20 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 40, J = N, 1, -1 IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) IX = JX DO 30, I = J - 1, 1, -1 IX = IX - INCX X( IX ) = X( IX ) - X( JX )*A( I, J ) 30 CONTINUE END IF JX = JX - INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = 1, N IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) DO 50, I = J + 1, N X( I ) = X( I ) - X( J )*A( I, J ) 50 CONTINUE END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) IX = JX DO 70, I = J + 1, N IX = IX + INCX X( IX ) = X( IX ) - X( JX )*A( I, J ) 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF END IF ELSE * * Form X := INV( A' )*X. * IF( LSAME( UPLO, 'U' ) )THEN IF( INCX.EQ.1 )THEN DO 100, J = 1, N DO 90, I = 1, J - 1 X( J ) = X( J ) - A( I, J )*X( I ) 90 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) 100 CONTINUE ELSE JX = KX DO 120, J = 1, N IX = KX DO 110, I = 1, J - 1 X( JX ) = X( JX ) - A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) JX = JX + INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = N, 1, -1 DO 130, I = N, J + 1, -1 X( J ) = X( J ) - A( I, J )*X( I ) 130 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )/A( J, J ) 140 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 160, J = N, 1, -1 IX = KX DO 150, I = N, J + 1, -1 X( JX ) = X( JX ) - A( I, J )*X( IX ) IX = IX - INCX 150 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) JX = JX - INCX 160 CONTINUE END IF END IF END IF * RETURN * * End of STRSV . * END LOGICAL FUNCTION LSAME ( CA, CB ) * .. Scalar Arguments .. CHARACTER*1 CA, CB * .. * * PURPOSE * ======= * * LSAME Tests if CA is the same letter as CB regardless of case. * CB is assumed to be an upper case letter. LSAME returns .TRUE. if * CA is either the same as CB or the equivalent lower case letter. * * 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. * * PARAMETERS * ========== * * 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 * 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 ) * .. Scalar Arguments .. INTEGER INFO CHARACTER*6 SRNAME * .. * * PURPOSE * ======= * * XERBLA Is an error handler for the level 2 blas routines. * * It is called by the level 2 blas routines if an input parameter is * invalid. * * Installers should 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. * * * AUXILIARY ROUTINE FOR LEVEL 2 BLAS. * * Written on 20-July-1986. * * .. Executable Statements .. * WRITE (6,99999) SRNAME, INFO * RETURN * 99999 FORMAT ( ' ** ON ENTRY TO ', A6, ' PARAMETER NUMBER ', I2, $ ' HAD AN ILLEGAL VALUE' ) * * End of XERBLA. * END End of hssxev.f echo index 1>&2 cat >index <<'End of index' From rgrimes@BOEING.COM Thu Oct 29 22:36:16 1987 Date: Thu, 29 Oct 87 08:52:40 pst From: Roger Grimes Subject: out-of-core eigensolver 3 files the certification program , the source code, and the tex documentation for the out-of-core eigensolver hssxev - An out-of-core symmetric eigensolve r for large dense problems. It uses block householder reductions to reduce the full dense matrix to banded form. The banded form is then reduced to tridiagonal form and all eigenvalues are computed. Specified eigenvectors are computed using inverse iteration with the band matrix and then back transformed to orginal form. End of index echo test-program.f 1>&2 cat >test-program.f <<'End of test-program.f' PROGRAM CSSXEV C C----------------------------------------------------------------------- C C CSSXEV IS THE MAIN PROGRAM TO TEST THE COMPUTATION OF C EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MATRIX C WHICH IS TOO LARGE TO FIT IN CENTRAL MEMORY. C C WRITTEN BY ROGER G. GRIMES, BOEING COMPUTER SERVICES, BELLEVUE, C WASHINGTON, JANUARY, 1987. C C----------------------------------------------------------------------- C IMPLICIT COMPLEX ( A - Z ) C C --------------------------- C ... PROBLEM SIZE PARAMETERS C --------------------------- C INTEGER BLKSIZ, NEQNS, NUMVEC, NWORK C PARAMETER ( BLKSIZ = 50, NEQNS = 500, NUMVEC = 10, 2 NWORK = 20 + 3 * NEQNS * BLKSIZ + NEQNS 3 + BLKSIZ + BLKSIZ**2 ) C C----------------------------------------------------------------------- C C ------------------------- C ... VARIABLE DECLARATIONS C ------------------------- C INTEGER I, IERR, J, L, OUTPUT, 1 UNIT1, UNIT2, UNIT3 C INTEGER INDEX(NUMVEC) C LOGICAL PASS C REAL MF, T, T1, T2 C REAL A(NEQNS,NEQNS), EIGVEC(NEQNS,NUMVEC), 1 EIGVAL(NEQNS), WORK(NWORK), 3 NUMER(NUMVEC,NUMVEC), DENOM(NUMVEC,NUMVEC) C REAL OPS COMMON / RGGOPS / OPS C C -------------------- C ... SUBPROGRAMS USED C -------------------- C REAL RANF, SDOT C EXTERNAL HSSXEV, HSSXE1, PRTVEC, RANF, 1 SDOT, SECOND, SGEMV C C ------------------- C ... DATA STATEMENTS C ------------------- C DATA OUTPUT / 6 /, 1 UNIT1 / 1 /, UNIT2 / 2 /, UNIT3 / 3 / C C----------------------------------------------------------------------- C WRITE ( OUTPUT, 20000 ) NEQNS, BLKSIZ, NUMVEC, NWORK 20000 FORMAT ( '1' ///5X, ' TESTING OF SUBROUTINE HSSXEV' 1 /10X, 'NUMBER OF EQUATIONS = ', I10 2 /10X, 'BLOCK SIZE = ', I10 4 /10X, 'NUMBER OF EIGENVECTORS = ', I10 5 /10X, 'SIZE OF WORKING STORAGE = ', I10 ) C IERR = 0 OPS = 0. C C ---------------------------- C ... GENERATE ORIGINAL MATRIX C ---------------------------- C DO 100 J = 1, NEQNS A(J,J) = RANF() DO 100 I = J+1, NEQNS A(I,J) = RANF() A(J,I) = A(I,J) 100 CONTINUE C C ------------------------------------- C ... USE HSSXEV TO COMPUTE EIGENVALUES C ------------------------------------- C OPS = 0. C CALL SECOND ( T1 ) C CALL HSSXEV ( 0, A, NEQNS, BLKSIZ, UNIT1, UNIT2, 1 WORK, NWORK, EIGVAL, IERR ) IF ( IERR .NE. 0 ) GO TO 8100 C DO 200 I = 1, NEQNS CALL HSSXEV ( 1, A(I,I), NEQNS, BLKSIZ, UNIT1, UNIT2, 1 WORK, NWORK, EIGVAL, IERR ) IF ( IERR .NE. 0 ) GO TO 8100 200 CONTINUE C CALL SECOND ( T2 ) T2 = T2 - T1 MF = OPS / 1.0E6 / T2 C WRITE ( OUTPUT, '(''0TIME AND RATE FOR HSSXEV = '', 2F12.6)') 1 T2, MF C C -------------------------- C ... PRINT OUT EIGENVALUES. C -------------------------- C L = MIN ( NEQNS, 50 ) C CALL PRTVEC ( 'EIGENVALUES COMPUTED USING HSSXEV', 1 L, EIGVAL, OUTPUT ) C C --------------------------------- C ... COMPUTE SELECTED EIGENVECTORS C --------------------------------- C DO 400 I = 1, NUMVEC INDEX(I) = I 400 CONTINUE C OPS = 0. C CALL SECOND ( T1 ) C CALL HSSXE1 ( NEQNS, BLKSIZ, UNIT1, UNIT2, UNIT3, EIGVAL, 1 NUMVEC, INDEX, WORK, NWORK, IERR ) C CALL SECOND ( T2 ) T2 = ( T2 - T1 ) / NUMVEC MF = OPS / 1.0E6 / NUMVEC / T2 C WRITE ( OUTPUT, '(''0TIME AND RATE FOR HSSXE1 PER VECTOR = '', 1 2F12.6)') T2, MF C IF ( IERR .NE. 0 ) GO TO 8200 C C ------------------------------------- C ... READ IN THE COMPUTED EIGENVECTORS C ------------------------------------- C DO 500 J = 1, NUMVEC READ ( UNIT3, REC=J ) ( EIGVEC(I,J), I = 1, NEQNS ) 500 CONTINUE C C -------------------------------------------------------- C ... COMPUTE THE MATRICES NUMER = X'*A*X AND DENOM = X'*X C -------------------------------------------------------- C DO 610 J = 1, NUMVEC C C CALL HSMVPS ( 1, NEQNS, NEQNS, A, NEQNS, C 1 A(1,J), WORK, IERR) CALL SGEMV ( 'N', NEQNS, NEQNS, 1.0, A, NEQNS, 1 EIGVEC(1,J), 1, 0.0, WORK, 1 ) C DO 600 I = 1, NUMVEC C NUMER(I,J) = SDOT ( NEQNS, EIGVEC(1,I), 1, 1 WORK, 1 ) DENOM(I,J) = SDOT ( NEQNS, EIGVEC(1,I), 1, 1 EIGVEC(1,J), 1 ) C 600 CONTINUE C 610 CONTINUE C C ------------------------------------- C ... CHECK THE RAYLEIGH-RITZ QUOTIENTS C ------------------------------------- C PASS = .TRUE. C DO 700 I = 1, NUMVEC C T = ABS ( EIGVAL(INDEX(I)) - NUMER(I,I) ) 1 / MAX ( ABS ( EIGVAL(INDEX(I)) ), 1.0E0 ) C IF ( T .GT. 1.0E-10 ) THEN PASS = .FALSE. WRITE ( OUTPUT, 2200 ) I, T 2200 FORMAT ( ' *** FATAL ERROR *** RAYLEIGH-RITZ ', 1 'TEST FAILED FOR I = ', I5, 2 ' A-INNER PRODUCT = ', 1PE12.3 ) END IF C 700 CONTINUE C C ------------------------------------------------ C ... CHECK THE ORTHONORMALITY OF THE EIGENVECTORS C ------------------------------------------------ C DO 810 J = 1, NUMVEC C DO 800 I = 1, NUMVEC C IF ( I .EQ. J ) THEN T = ABS ( 1.0E0 - DENOM(I,J) ) ELSE T = ABS ( DENOM(I,J) ) END IF C IF ( T .GT. 1.0E-10 ) THEN PASS = .FALSE. WRITE ( OUTPUT, 2300 ) I, J, T 2300 FORMAT ( ' *** FATAL ERROR *** EIGENVECTOR ', 1 'ORTHONORMALITY TEST FAILED FOR I, J, = ', 2 2I5, ' INNER PRODUCT = ', 1PE12.3 ) END IF C 800 CONTINUE C 810 CONTINUE C IF ( PASS ) THEN WRITE ( OUTPUT, 21000 ) 21000 FORMAT ( ' HSSXEV ** PASSED ** CERTIFICATION' ) ELSE WRITE ( OUTPUT, 22000 ) 22000 FORMAT ( ' HSSXEV ** FAILED ** CERTIFICATION' ) END IF C GO TO 9000 C C----------------------------------------------------------------------- C C --------------- C ... ERROR TRAPS C --------------- C C -------------------------------- C ... ERROR FROM SUBROUTINE HSSXEV C -------------------------------- C 8100 WRITE ( OUTPUT, 28100 ) IERR 28100 FORMAT ( 1X, '*** FATAL ERROR *** UNEXPECTED ERROR RETURN FROM ', 1 'SUBROUTINE HSSXEV - IERR = ', I5 ) GO TO 9000 C C -------------------------------- C ... ERROR FROM SUBROUTINE HSSXE1 C -------------------------------- C 8200 WRITE ( OUTPUT, 28200 ) IERR 28200 FORMAT ( 1X, '*** FATAL ERROR *** UNEXPECTED ERROR RETURN FROM ', 1 'SUBROUTINE HSSXE1 - IERR = ', I5 ) GO TO 9000 C C----------------------------------------------------------------------- C C ------------------------- C ... END OF PROGRAM CSSXEV C ------------------------- C 9000 CONTINUE STOP END SUBROUTINE PRTVEC ( TITLE, N, X, OUTPUT ) C C----------------------------------------------------------------------- C C PRTVEC PRINTS OUT THE VECTOR X OF LENGTH N TO OUTPUT LOGICAL C UNIT OUTPUT USING THE CHARACTER STRING IN TITLE AS A TITLE. C C WRITTEN BY ROGER G. GRIMES, BOEING COMPUTER SERVICES, BELLEVUE, C WASHINGTON, JANUARY, 1987. C C----------------------------------------------------------------------- C C ------------------------ C ... VARIABLE DECLARATION C ------------------------ C IMPLICIT COMPLEX ( A - Z ) C INTEGER I, L, N, OUTPUT C REAL X(N) C CHARACTER*(*) TITLE CHARACTER* 80 LINE C C----------------------------------------------------------------------- C C --------------- C ... WRITE TITLE C --------------- C L = MIN ( LEN ( TITLE ), 80 ) C DO 100 I = 1, L LINE(I:I) = '-' 100 CONTINUE C DO 200 I = L+1, 80 LINE(I:I) = ' ' 200 CONTINUE C WRITE ( OUTPUT, 2000 ) TITLE, LINE(1:L) 2000 FORMAT ( /1X, A /1X, A / ) C C ------------------------ C ... WRITE OUT THE VECTOR C ------------------------ C WRITE ( OUTPUT, 2100 ) X 2100 FORMAT ( 5X, 1P5E15.6 ) C C----------------------------------------------------------------------- C C ----------------- C ... END OF PRTVEC C ----------------- C RETURN END End of test-program.f echo tex 1>&2 cat >tex <<'End of tex' \newpage \subprogramtitle{HSSXEV, HSSXE1 - Eigenvalues/Eigenvectors for Real Symmetric Eigenproblems with External Storage} \versions HSSXEV, HSSXE1 - REAL \endversions \purpose HSSXEV finds all of the eigenvalues of a REAL symmetric, dense eigenproblem $AX = X \Lambda$. HSSXE1 is then used to compute the eigenvectors for selected eigenvalues. HSSXEV and HSSXE1 are intended for use for matrices which are so large as to require external (``out-of-core'') storage for $A$. BCSLIB Subroutine HSSYEV is generally recommended for problems which are small enough for internal solution. \endpurpose \method HSSXEV first uses block Householder transformations to reduce $A$ to a band matrix which can be stored internally and whose eigenvalues are identical to the eigenvalues of $A$. Then Givens plane rotations are used to reduce the banded matrix to tridiagonal form while maintaining the similarity of the eigenvalues. Finally the eigenvalues of the tridiagonal matrix, which are also the eigenvalues of $A$, are computed using the $QL$ algorithm. HSSXE1 computes the eigenvectors of $A$ for user specified eigenvalues by first computing the eigenvectors of the band matrix using inverse iteration. The eigenvectors of the band matrix are then back transformed to the eigenvectors of $A$. The user calls HSSXEV repetitively to compute the eigenvalues of $A$. The first call is for initialization which initializes the working array and opens I/O files on external storage. The user must specify both the matrix order $n$ and the block size $p$. The next $n$ calls pass the rows of the upper triangular portion of $A$ to HSSXEV. HSSXEV collects the matrix $A$ into its working array and when $p$ rows are acculumated they are written to external storage. During the call which inputs the last row of $A$, the matrix is reduced and the eigenvalues are computed and returned. At the completion of HSSXEV, the band matrix and the block Householder transformations are on external storage. Since $A$ is symmetric, $A$ can be input by columns of the lower triangle as well. For more information, see REMARKS and EXAMPLE. After the last call to HSSXEV and the eigenvalues have been computed, HSSXE1 can be called to compute eigenvectors corresponding to selected eigenvalues. The eigenvectors are written to external storage. HSSXE1 can be called several times to compute the eigenvectors corresponding to several sets of eigenvalues. For more information, see REMARKS and EXAMPLE. The processing by HSSXEV is directed by MODE, an initialization/input flag. The user sets MODE = 0 for the initialization call and MODE = 1 for each of the $n$ input calls. \endmethod \subprogramtitle{HSSXEV - Eigenvalues for Real Symmetric Matrices with External Storage} \purpose HSSXEV finds all of the eigenvalues of a REAL symmetric, dense eigenproblem $AX = X \Lambda$ using external storage. \endpurpose \usage REAL A(N),HOLD(NHOLD) CALL HSSXEV (MODE,A,N,P,LUNIT1,LUNIT2,HOLD,NHOLD,EIGVAL,IER) \endusage \beginarguments{MODE} Initialization/input flag, see METHOD for explanation. \splitline{MODE=0} Initialization. $A$ not used. \splitline{MODE=1} $A$ is an input argument from which HSSXEV obtains a row of the upper triangular portion of the matrix. \inp{A} If MODE=1, A is an array of length $N-I+1$ which contains the upper triangular part of row $I$ of the matrix to be input. The diagonal entry of the matrix is stored in position 1 of A while the entry for the last column of the matrix is stored in position $N-I+1$. It is used only when MODE=1. EXAMPLE illustrates how to call HSSXEV so that successive A arrays correspond to successive rows in the upper triangular portion of the matrix to be factored. \inp{N} Order of the matrix $A$, N $>$ 0. \inp{P} Block size, P $>$ 0. \inp{LUNIT1} Logical unit number of an input/output dataset which will be opened and used as a direct-access dataset to store the matrix $A$ and is overwritten with the band matrix. The approximate total amount of I/O transfer to this dataset will be $n^3/p$ real numbers. The length of the dataset will be approximately $n^2/2$ real numbers. See REMARKS. \inp{LUNIT2} Logical unit number of an input/output dataset which will be opened and used as a direct-access dataset by HSSXEV to store the block Householder transformations for later use by HSSXE1. The approximate total amount of I/O transfer by HSSXEV to this dataset will be $n^2$ real numbers. The approximate total amount of I/O transfer for each call to HSSXE1 from this dataset will be $n^2$ real numbers. The length of the dataset will be approximately $n^2$ real numbers. See REMARKS. \beginworkingstorage{HOLD} Work array of length NHOLD. This array must be preserved from the initial call (MODE=0) to HSSXEV until all desired eigenvectors have been computed. \workingstorage{NHOLD} The length of array HOLD. NHOLD must be at least $$ 20 + 3\rm{NP} + \rm{N} + \rm{P} + max ( \rm{P}^2, 2\rm{N} ). $$ \beginoutput{EIGVAL} If IER $\ge$ 0, array EIGVAL contains the computed eigenvalues in ascending order. If IER = 0 all eigenvalues have been computed successfully. If IER $>$ 0 then only the first IER-1 eigenvalues have been computed. If IER $<$ 0, EIGVAL is unchanged from input. \out{IER} Success/error code. Results have not been computed for IER $<$ 0. \splitline{IER=0} Success, all eigenvalues have been computed. \splitline{IER$>$0} Only the first IER-1 eigenvalues have been computed successfully. \splitline{IER=$-$1} Invalid value for MODE. \splitline{IER=$-$2} N $\le$ 0. \splitline{IER=$-$3} P $\le$ 0 or P $>$ N/4. \splitline{IER=$-$4} NHOLD too small. \splitline{IER=$-$5} I/O error on LUNIT1. \splitline{IER=$-$6} I/O error on LUNIT2. \splitline{IER=$-$7} Unexpected error return from subroutine \break HSMMPG. \endarguments \subprogramtitle{HSSXE1 - Compute Selected Eigenvectors for Real Symmetric Matrices with External Storage} \purpose HSSXE1 computes eigenvectors for specified eigenvalues after HSSXEV has been used to find all of the eigenvalues of a REAL general, dense matrix $A$ using external storage. \endpurpose \usage REAL A(N),HOLD(NHOLD) CALL HSSXEV (MODE,A,N,P,LUNIT1,LUNIT2,HOLD,NHOLD,EIGVAL,IER) CALL HSSXE1 (N,P,LUNIT1,LUNIT2,LUNIT3,EIGVAL,NUMVEC,INDEX, \continue{1}{\blank{13}HOLD,NHOLD,IER)} \endusage \beginarguments{N} Order of the matrix $A$, N $>$ 0. \inp{P} Block size, P $>$ 0. \inp{LUNIT1} Logical unit number of the first input/output dataset supplied to HSSXEV. \inp{LUNIT2} Logical unit number of the second input/output dataset supplied to HSSXEV. \inp{LUNIT3} Logical unit number of an output dataset which will be opened and used as a direct-access dataset by HSSXE1 to store the computed eigenvectors. The total amount of I/O transfer by HSSXE1 to this dataset will be $n \times$ NUMVEC real numbers. The length of the dataset is also $n \times$ NUMVEC real numbers. See REMARKS. \inp{EIGVAL} Array containing the eigenvalues computed by HSSXEV. \inp{NUMVEC} Number of eigenvectors to compute, NUMVEC $>$ 0. \inp{INDEX} Array containing the indices of the eigenvalues in EIGVAL for which the associated eigenvectors are to be computed. \beginworkingstorage{HOLD} Work array of length NHOLD. This array must be preserved from the initial call (MODE=0) to HSSXEV until all desired eigenvectors have been computed. \workingstorage{NHOLD} The length of array HOLD. NHOLD must be at least $$ 20 + 3\rm{NP} + \rm{N} + \rm{P} + max ( \rm{P}^2, 2\rm{N} ). $$ \beginoutput{IER} Success/error code. Results have not been computed for IER $<$ 0. \splitline{IER=0} Success, all eigenvectors have been computed. \splitline{IER$>$0} Only the first IER-1 eigenvectors have been computed successfully. \splitline{IER=$-$1} Illegal processing sequence. \splitline{IER=$-$2} N $\le$ 0. \splitline{IER=$-$3} P $\le$ 0 or P $>$ N/4. \splitline{IER=$-$4} NUMVEC $\le$ 0. \splitline{IER=$-$5} NHOLD too small. \splitline{IER=$-$6} I/O error on LUNIT1. \splitline{IER=$-$7} I/O error on LUNIT2. \splitline{IER=$-$8} I/O error on LUNIT3. \splitline{IER=$-$9} Unexpected error return from subroutine \break HSMMPG. \endarguments \remarks For each problem, the same values of N, P, LUNIT1, and LUNIT2 and the same HOLD array must be used in all calls to HSSXEV. When computing eigenvectors with HSSXE1 the values of N, P, LUNIT1, and LUNIT2 as well as the same arrays for HOLD and EIGVAL must be used as in the call to HSSXEV. LUNIT3, NUMVEC, and INDEX may vary in calling HSSXE1. HSSXEV uses a block algorithm where the block size, $p$, is provided by the user to HSSXEV. The choice of block size is dependent on the matrix order, $n$, and affects the length of the HOLD array, NHOLD. The number of words of actual workspace used, given $n$ and $p$, is $$ \hfil 20 + 3np + n + p + max ( p^2, 2n ). \hfil $$ The smallest blocksize allowed is $p = 1$. However the efficiency of the algorithm is severly hindered by block sizes which are too small or too large. A block size of 50 to 60 is recommended for obtain near optimal efficiency. This requires NHOLD to be from $153n$ to $183n$. The algorithm requires $p \le n/4$. Logical unit number LUNIT1 is first used to store the matrix $A$. This matrix is successively transformed until it becomes a banded matrix. At the successful completion of HSSXEV the banded matrix is stored on LUNIT1. The amount of I/O transfer to LUNIT1 is approximately $n^3/p$ real words while the length of the dataset is $n^2/2$ real words. Logical unit number LUNIT2 is use to store the block Householder transformations required to transform $A$ to banded form. These transformations are used by HSSXE1 to back transform the eigenvectors of the banded matrix to those of $A$. The amount of I/O transfer to LUNIT2 is approximately $n^2$ real words at the completion of HSSXEV and $n^2$ additional words for each set of $p$ eigenvectors computed in HSSXE1. The length of the dataset is $n^2$ real words. Logical unit number LUNIT3 is used to store the eigenvectors computed in a single call to HSSXE1. It is a FORTRAN direct access file with each record having record length equivalent to $n$ real words. The eigenvectors are written in the same order as the specification in the INDEX array. The eigenvectors can be read from LUNIT3 after completion of a successful HSSXE1 using the standard FORTRAN-77 direct access READ statement. The $i$-th record contains the eigenvectors which corresponds to the eigenvalue specified by INDEX($i$). See EXAMPLE. On some systems, it may be necessary to allocate a sufficient amount of space to LUNIT1, LUNIT2, and LUNIT3 before the initialization call for a large problem. It may also be desirable to assign the I/O units to a particular type of I/O device. The I/O units are not closed at the end of HSSXEV or HSSXE1. If the user wishes to reuse them for another problem he must close them and then dispose of them before reusing them. The user should not manipulate LUNIT1 and LUNIT2 while HSSXEV is active or while additional eigenvectors are left to be computed. \endremarks \vfill \pagebreak \titledremarks{EXAMPLE}Compute the eigenvalues for a REAL symmetric eigenproblem of order 1000. Compute the eigenvectors corresponding to all eigenvalues in the interval {\lbrack -10., 10.\rbrack}. NHOLD is large enough for a block size of 50. \code INTEGER INDEX(1000) REAL A(1000),EIGVAL(1000),HOLD(153570),EIGVEC(1000,100) DATA NHOLD/153570/ N = 1000 P = 50 LUNIT1 = 1 LUNIT2 = 2 LUNIT3 = 3 C C ... INITIALIZATION CALL TO HSSXEV C MODE = 0 CALL HSSXEV (MODE,A,N,P,LUNIT1,LUNIT2,HOLD,NHOLD,EIGVAL,IER) IF (IER.NE.0) GO TO 100 C C ... INPUT MATRIX TO HSSXEV. LAST CALL WILL COMPUTE THE C EIGENVALUES. C MODE = 1 DO 10 I = 1, N (Evaluate row I of $A$. Store the entries on and to the left of the main diagonal in positions 1 through N-I+1 of A.) CALL HSSXEV (MODE,A,N,P,LUNIT1,LUNIT2,HOLD,NHOLD,EIGVAL, 1 IER) IF (IER.NE.0) GO TO 100 10 CONTINUE C C ... SELECT THE EIGENVALUES WHOSE EIGENVECTORS ARE C TO BE COMPUTED C NUMVEC = 0 DO 20 I = 1, N IF ( EIGVAL(I) .GE. -10. .AND. EIGVAL(I) .LE. 10. ) THEN NUMVEC = NUMVEC + 1 INDEX(NUMVEC) = I END IF 20 CONTINUE \vfill \pagebreak C C ... COMPUTE THE EIGENVECTORS C CALL HSSXE1 (N,P,LUNIT1,LUNIT2,LUNIT3,HOLD,NHOLD,EIGVAL, 1 NUMVEC,INDEX,IER) IF (IER.NE.0) GO TO 100 C C ... READ IN THE COMPUTED EIGENVECTORS C DO 30 I = 1, NUMVEC READ ( LUNIT3, REC=I ) ( EIGVEC(J,I), J = 1, N ) 30 CONTINUE . {\rm (Continue processing)} . 100 {\rm (Error processing)} \endcode \endtitledremarks End of tex