*********************************************************************** * The FORTRAN source code for the generic version of LANZ. * * * * The first few routines in this file comprise the test program. * * They are separated from the rest of the file by a comment * * section similar to this one. The rest of the file contains * * the LANZ source code. * * * * For illustration purposes, the installation steps needed to * * create an executable verion of LANZ, including the test * * programs, on the IBM 3090 running AIX are now given: * * First the two files are given suffixes that reflect the type * * of source code contained in them: * * mv lanz.file3 lanz.file3.f * * mv lanz.file4 lanz.file4.c * * * * Next the C compiler is invoked to create lanz.file4.o: * * cc -O -c lanz.file4.c * * * * Lastly the FORTRAN compiler is called to create the * * executable, lanz.ibm: * * fvs -xa -f'opt(3) vec' -o lanz.ibm lanz.file3.f lanz.file4.o * * * *********************************************************************** C$FORTRAN BISEC1 *********************************************************************** * LANZ SOFTWARE PACKAGE * * FILENAME: BISECT.MSC * * AUTHOR: MARK JONES * * PURPOSE: USED IN THE TRIDIAGONAL EIGENSOLVER TO PROTECT THE * * NEWTON ROOTFINDER FROM STRAYING OUT OF BOUNDS * *********************************************************************** SUBROUTINE BISEC1(J,ALF,BET2,LEFT,RIGHT,EPS,INDEX,DEBUG) * THE ORDER OF T(J) INTEGER J * ALF IS THE DIAGONAL OF T(J) DOUBLE PRECISION ALF(J) * BET2 IS THE SQUARE OF THE OFF-DIAGONAL OF T(J) DOUBLE PRECISION BET2(J) * LEFT (RIGHT) BOUND ON INTERVAL DOUBLE PRECISION LEFT, RIGHT * EPSILON IS THE MACHINE PRECISION DOUBLE PRECISION EPS * THE NUMBER OF THE EIGENVALUE THAT WE ARE SEEKING INTEGER INDEX * THE LEVEL OF DEBUGGING OUTPUT INTEGER DEBUG * THE FOLLOWING ARE INTERNAL VARIABLES * THE SIZE OF THE BOUND DOUBLE PRECISION DIST, ODIST * THE DIRECTION CURRENTLY MOVING INTEGER DIRECT * HAVE WE SWITCHED DIRECTION? LOGICAL BTRACK * NUMBERS OF EIGENVALUES AT CERTAIN PROBES INTEGER NUMRIG, NUMLEF, NUMMID * THE MID POINT OF THE RANGE DOUBLE PRECISION MID * AN INTEGER FUNCTION INTEGER NUMLES EXTERNAL NUMLES ODIST = RIGHT - LEFT IF (ODIST.LE.0.0D0) THEN ODIST = ABS(RIGHT * (4.0D0*EPS)) IF (ODIST.EQ.0.0D0) THEN ODIST = EPS ENDIF ENDIF DIST = ODIST DIRECT = 0 BTRACK = .FALSE. * ADJUST THE RIGHT SIDE OF THE BOUND 10 CONTINUE NUMRIG = NUMLES(ALF,BET2,RIGHT,J,1,EPS) IF (NUMRIG.LT.INDEX) THEN IF (DEBUG.GT.0) THEN PRINT *,'WARNING: BISECT DETECTS NUMERICAL PROBLEM' PRINT *,'FIXED:ADJUSTED RIGHT SIDE TO RIGHT' PRINT *,DIRECT,' ',DIST,' ',ODIST,' ',RIGHT,' ',J ENDIF RIGHT = RIGHT + DIST IF ((DIRECT.EQ.2).OR.(BTRACK)) THEN DIST = DIST/2 IF (DIST.LE.0.0D0) THEN DIST = ABS(RIGHT * (2.0D0*EPS)) ENDIF BTRACK = .TRUE. ELSE DIST = DIST*2 ENDIF DIRECT = 1 GOTO 10 ELSE IF (NUMRIG.GT.INDEX) THEN IF (DEBUG.GT.0) THEN PRINT *,'WARNING: BISECT DETECTS NUMERICAL PROBLEM' PRINT *,'FIXED:ADJUSTED RIGHT SIDE TO LEFT' ENDIF RIGHT = RIGHT - DIST IF ((DIRECT.EQ.1).OR.(BTRACK)) THEN DIST = DIST/2 IF (DIST.LE.0.0D0) THEN DIST = ABS(RIGHT * (2.0D0*EPS)) ENDIF BTRACK = .TRUE. ELSE DIST = DIST*2 ENDIF DIRECT = 2 GOTO 10 ENDIF * ADJUST THE LEFT SIDE OF THE BOUND DIST = ODIST DIRECT = 0 BTRACK = .FALSE. 20 CONTINUE NUMLEF = NUMLES(ALF,BET2,LEFT,J,1,EPS) IF (NUMLEF.LT.INDEX-1) THEN IF (DEBUG.GT.0) THEN PRINT *,'WARNING: BISECT DETECTS NUMERICAL PROBLEM' PRINT *,'FIXED:ADJUSTED LEFT SIDE TO RIGHT' ENDIF LEFT = LEFT + DIST IF ((DIRECT.EQ.2).OR.(BTRACK)) THEN DIST = DIST/2 IF (DIST.LE.0.0D0) THEN DIST = ABS(LEFT * (2.0D0*EPS)) ENDIF BTRACK = .TRUE. ELSE DIST = DIST*2 ENDIF DIRECT = 1 GOTO 20 ELSE IF (NUMLEF.GT.INDEX-1) THEN IF (DEBUG.GT.0) THEN PRINT *,'WARNING: BISECT DETECTS NUMERICAL PROBLEM' PRINT *,'FIXED:ADJUSTED LEFT SIDE TO LEFT' ENDIF LEFT = LEFT - DIST IF ((DIRECT.EQ.1).OR.(BTRACK)) THEN DIST = DIST/2 IF (DIST.LE.0.0D0) THEN DIST = ABS(LEFT * (2.0D0*EPS)) ENDIF BTRACK = .TRUE. ELSE DIST = DIST*2 ENDIF DIRECT = 2 GOTO 20 ENDIF MID = (RIGHT + LEFT)/2.0D0 NUMMID = NUMLES(ALF,BET2,MID,J,1,EPS) IF (NUMMID.EQ.NUMRIG) THEN RIGHT = MID ELSE LEFT = MID ENDIF RETURN END C$FORTRAN BSTATE *********************************************************************** * LANZ SOFTWARE PACKAGE * * FILENAME: BOUND.MSC * * AUTHOR: MARK JONES * * PURPOSE: USED IN MSHIFT() TO DETERMINE THE BOUNDARIES WHEN * * SEARCHING FOR EIGENPAIRS IN A PRE-DETERMINED RANGE * *********************************************************************** SUBROUTINE BSTATE(STATE,ISLFT,ISRGT,LFTVAL,RGTVAL,RETVAL, C BLEFT,BRIGHT,RUN) * THE CURRENT STATE INTEGER STATE * DO WE HAVE LEFT AND RIGHT SHIFTS LOGICAL ISRGT, ISLFT * LEFT AND RIGHT SHIFTS DOUBLE PRECISION RGTVAL, LFTVAL * RETURN FROM LAST LANCZOS RUN INTEGER RETVAL * LEFT AND RIGHT BOUNDARIES DOUBLE PRECISION BLEFT,BRIGHT * THE LOOP VARIABLE FROM THE CALLING ROUTINE TO RESET INTEGER RUN IF (STATE.EQ.0) THEN RETURN ELSE IF (STATE.EQ.-1) THEN STATE = 1 ELSE IF (STATE.EQ.1) THEN IF (RETVAL.EQ.0) THEN STATE = 4 ISRGT = .TRUE. ISLFT = .TRUE. RGTVAL = BRIGHT LFTVAL = BLEFT RUN = 1 ELSE STATE = 2 ENDIF ELSE IF ((STATE.EQ.2).OR.(STATE.EQ.3)) THEN IF ((RGTVAL.GE.BRIGHT).OR.(RGTVAL.LE.BLEFT)) THEN ISRGT = .FALSE. ENDIF IF ((LFTVAL.GE.BRIGHT).OR.(LFTVAL.LE.BLEFT)) THEN ISLFT = .FALSE. ENDIF IF (.NOT.((ISRGT).OR.(ISLFT))) THEN STATE = 4 ISRGT = .TRUE. ISLFT = .TRUE. RGTVAL = BRIGHT LFTVAL = BLEFT ENDIF ELSE IF (STATE.EQ.4) THEN STATE = 5 ELSE IF (STATE.EQ.5) THEN STATE = 0 ENDIF RETURN END INTEGER FUNCTION GBOUND(INTIDX,SIGLST,INTLST,BLEFT,BRIGHT,JOB) * THE NUMBER OF INERTIA CHECKS INTEGER INTIDX * THE SIGMAS DOUBLE PRECISION SIGLST(*) * THE INERTIAS INTEGER INTLST(*) * THE LEFT AND RIGHT BOUNDS DOUBLE PRECISION BLEFT, BRIGHT * THE JOB TO DO INTEGER JOB * INTERNAL VARIABLES INTEGER NLEFT, NRIGHT, NMID DOUBLE PRECISION MID INTEGER I MID = (BRIGHT + BLEFT) / 2.0D0 DO 10 I = 1, INTIDX IF (SIGLST(I).EQ.BLEFT) THEN NLEFT = INTLST(I) ELSE IF (SIGLST(I).EQ.BRIGHT) THEN NRIGHT = INTLST(I) ELSE IF (SIGLST(I).EQ.MID) THEN NMID = INTLST(I) ENDIF 10 CONTINUE IF (JOB.EQ.0) THEN GBOUND = NRIGHT - NLEFT ELSE IF (JOB.EQ.4) THEN GBOUND = NRIGHT - NMID ELSE GBOUND = NMID - NLEFT ENDIF RETURN END C$FORTRAN BSORT *********************************************************************** * LANZ SOFTWARE PACKAGE * * FILENAME: BSORT.MSC * * AUTHOR: MARK JONES * * PURPOSE: BUBBLESORT ROUTINES USED TO SORT SMALL VECTORS * *********************************************************************** * SORTS 1 DOUBLE PRECISION VECTOR SUBROUTINE BSORT(N,INVEC) * THE LENGTH OF THE VECTOR INTEGER N * THE DOUBLE PRECISION VECTOR TO BE SORTED DOUBLE PRECISION INVEC(*) * THE REST ARE INTERNAL VARIABLES * COUNT VARIABLES INTEGER I,J * A TEMPORARY VARIABLE DOUBLE PRECISION TEMP DO 10 I = 1, N-1 DO 20 J = I+1, N IF (INVEC(I).GT.INVEC(J)) THEN TEMP = INVEC(I) INVEC(I) = INVEC(J) INVEC(J) = TEMP ENDIF 20 CONTINUE 10 CONTINUE RETURN END * SORTS A DOUBLE PRECISION VECTOR WITH A DATA INTEGER VECTOR SUBROUTINE BSORT2(N,KEYVEC,DVEC) * THE LENGTH OF THE VECTOR INTEGER N * THE DOUBLE PRECISION VECTOR TO BE SORTED DOUBLE PRECISION KEYVEC(*) * THE INTEGER VECTOR ASSOCIATED WITH KEYVEC INTEGER DVEC(*) * THE REST ARE INTERNAL VARIABLES * COUNT VARIABLES INTEGER I,J * A TEMPORARY VARIABLE DOUBLE PRECISION TEMP * A TEMPORARY VARIABLE INTEGER ITEMP DO 10 I = 1, N-1 DO 20 J = I+1, N IF (KEYVEC(I).GT.KEYVEC(J)) THEN TEMP = KEYVEC(I) KEYVEC(I) = KEYVEC(J) KEYVEC(J) = TEMP ITEMP = DVEC(I) DVEC(I) = DVEC(J) DVEC(J) = ITEMP ENDIF 20 CONTINUE 10 CONTINUE RETURN END * SORTS A DOUBLE PRECISION VECTOR WITH TWO DATA VECTORS SUBROUTINE BSORT3(N,KEYVEC,DVEC1,DVEC2) * THE LENGTH OF THE VECTOR INTEGER N * THE DOUBLE PRECISION VECTOR TO BE SORTED DOUBLE PRECISION KEYVEC(*) * THE INTEGER VECTOR ASSOCIATED WITH KEYVEC INTEGER DVEC1(*) * THE DOUBLE PRECISION VECTOR ASSOCIATED WITH KEYVEC DOUBLE PRECISION DVEC2(*) * THE REST ARE INTERNAL VARIABLES * COUNT VARIABLES INTEGER I,J * A TEMPORARY VARIABLE DOUBLE PRECISION TEMP * A TEMPORARY VARIABLE INTEGER ITEMP DO 10 I = 1, N-1 DO 20 J = I+1, N IF (KEYVEC(I).GT.KEYVEC(J)) THEN TEMP = KEYVEC(I) KEYVEC(I) = KEYVEC(J) KEYVEC(J) = TEMP ITEMP = DVEC1(I) DVEC1(I) = DVEC1(J) DVEC1(J) = ITEMP TEMP = DVEC2(I) DVEC2(I) = DVEC2(J) DVEC2(J) = TEMP ENDIF 20 CONTINUE 10 CONTINUE RETURN END * SORTS AN INTEGER VECTOR WITH A DATA DOUBLE PRECISION VECTOR SUBROUTINE BSORTT(N,KEYVEC,DVEC) * THE LENGTH OF THE VECTOR INTEGER N * THE DOUBLE PRECISION VECTOR TO BE SORTED INTEGER KEYVEC(*) * THE INTEGER VECTOR ASSOCIATED WITH KEYVEC DOUBLE PRECISION DVEC(*) * THE REST ARE INTERNAL VARIABLES * COUNT VARIABLES INTEGER I,J * A TEMPORARY VARIABLE INTEGER TEMP * A TEMPORARY VARIABLE DOUBLE PRECISION ITEMP DO 10 I = 1, N-1 DO 20 J = I+1, N IF (KEYVEC(I).GT.KEYVEC(J)) THEN TEMP = KEYVEC(I) KEYVEC(I) = KEYVEC(J) KEYVEC(J) = TEMP ITEMP = DVEC(I) DVEC(I) = DVEC(J) DVEC(J) = ITEMP ENDIF 20 CONTINUE 10 CONTINUE RETURN END C$FORTRAN BUNCH *********************************************************************** * LANZ SOFTWARE PACKAGE * * FILENAME: BUNCH.MSC * * AUTHOR: MARK JONES * * PURPOSE: BANDED BUNCH-KAUFMAN FACTORIZATION AND FORWARD/BACK * * SOLUTION AS DESCRIBED BY JONES AND PATRICK * * CONTAINS BOTH A VECTOR AND A NON-VECTOR VERSION * *********************************************************************** * THE NON-VECTOR VERSION SUBROUTINE BUNCH(N,K,KRP,KCP,M,MRP,MCP, C A,ROWPTR,RWSTRT,B,X,SIGMA,INFLAG,MXSIZE, C DIAG,SUBDAG,V,V2,V3,V4,V5,V6,VL,VL2,VL3,VL4,VL5,VL6, C P,IUPTO,INERT,DEBUG,CMOFF,CMRP,CMCP, C CKOFF,CKRP,CKCP,UNROLL) * N IS THE ORDER OF THE MATRIX INTEGER N * K IS A SPARSE MATRIX DOUBLE PRECISION K(*) * KRP CONTAINS THE ROW POINTERS FOR K INTEGER KRP(*) * KCP CONTAINS THE COLUMN NUMBERS FOR K INTEGER KCP(*) * M IS A SPARSE MATRIX DOUBLE PRECISION M(*) * MRP CONTAINS THE ROW POINTERS FOR M INTEGER MRP(*) * MCP CONTAINS THE COLUMN NUMBERS FOR M INTEGER MCP(*) * A IS THE FACTORED VERSION OF (K-SIGMA*M) DOUBLE PRECISION A(*) * ROWPTR POINTS TO THE ROWS OF A INTEGER ROWPTR(*) * RWSTRT CONTAINS THE COLUMN NUMBER OF (K-SIGMA*M) WHERE * EACH ROW STARTS INTEGER RWSTRT(*) * B IS THE RIGHT-HAND SIDE DOUBLE PRECISION B(*) * X IS THE SOLUTION TO AX=B DOUBLE PRECISION X(*) * SIGMA IS THE SHIFT BEING USED DOUBLE PRECISION SIGMA * INDICATES WHAT COMBINATION OF FACT, FORWARD AND BACK SOLVE TO DO INTEGER INFLAG * THE MAXIMUM AMOUNT OF SPACE ALLOCATED FOR A INTEGER MXSIZE * WORK VECTORS FOR THE ALGORITHM * THE DIAGONAL OF THE D MATRIX DOUBLE PRECISION DIAG(*) * THE OFF-DIAGONAL OF THE D MATRIX DOUBLE PRECISION SUBDAG(*) * THE FOLLOWING VECTORS ARE USED WHEN UPDATING THE A MATRIX DOUBLE PRECISION V(*) DOUBLE PRECISION V2(*) DOUBLE PRECISION V3(*) DOUBLE PRECISION V4(*) DOUBLE PRECISION V5(*) DOUBLE PRECISION V6(*) DOUBLE PRECISION VL(*) DOUBLE PRECISION VL2(*) DOUBLE PRECISION VL3(*) DOUBLE PRECISION VL4(*) DOUBLE PRECISION VL5(*) DOUBLE PRECISION VL6(*) * A VECTOR TO KEEP TRACK OF PERMUTATIONS INTEGER P(*) * RECORDS HOW MANY COLUMNS CONTAIN NON-ZEROES BELOW A ROW INTEGER IUPTO(*) * THE INERT INTEGER INERT * THE DEBUGGING LEVEL INTEGER DEBUG * THE FOLLOWING ARE THE K AND M MATRICES STORED BY COLUMNS DOUBLE PRECISION CMOFF(*), CKOFF(*) INTEGER CMRP(*), CKRP(*) INTEGER CMCP(*), CKCP(*) * THE LEVEL OF LOOP UNROLLING INTEGER UNROLL * SOME INTERNAL VARIABLES * COUNT VARIABLES INTEGER L, I, J, START, RSTART, NSTART, ISTART, MYNUM, IMULT INTEGER IT1, IT2, IT3, IT4, IT5, IT6, IT7, IT8, II, IOFF * THE MAXIMUM LEVEL OF COMBINING OF ROWS IN THE TRIANGULAR SOLVES INTEGER ICMAX * THE NUMBER OF ROWS TO WORK ON IN THE TRIANGULAR SOLVES INTEGER NUMROW * RECORDS THE LAST ROW READ FROM K AND M INTO A INTEGER UPTO * THE ROW AT STEP L CONTAINING THE LARGEST ELEMENT * JUST AS IN THE JONES-PATRICK PAPER INTEGER TR(1), R * THE FIRST EMPTY ELEMENT IN A INTEGER NEXT * USED TO STORE THE RESULTS OF MIN AND MAX CALLS INTEGER TMPMAX * TEMPORARY VARIABLES FOR LOOPS DOUBLE PRECISION T, T1, T2, T3, T4, T5, T6, T7, TX * TEMPORARY VARIABLES FOR LOOPS DOUBLE PRECISION AT, AT1, AT2, AT3, AT4, AT5 * SOME OF THE VALUES IN THE PIVOT BLOCK * A * B C * D E F * G H I J * K L M N O * P Q R S T U DOUBLE PRECISION PAVAL, PBVAL, PCVAL, PDVAL, PEVAL, PFVAL, C PGVAL, PHVAL, PIVAL, PJVAL, PKVAL, PLVAL, PMVAL, PNVAL, POVAL, C PPVAL, PQVAL, PRVAL, PSVAL, PTVAL, PUVAL * TEMPORARY HOLDING VALUES FOR THE PIVOT BLOCK DOUBLE PRECISION TBVAL, TDVAL, TEVAL, TGVAL, THVAL, TIVAL, C TKVAL, TLVAL, TMVAL, TNVAL, TPVAL, TQVAL, TRVAL, TSVAL, TTVAL * STABILITY TEST LEVELS FOR PIVOTS DOUBLE PRECISION SLVL1, SLVL2, SLVL3, SLVL4, SLVL5, SLVL6, SLVLM * THE ELEMENT GROWTH BOUNDS FOR PIVOTS DOUBLE PRECISION GRATE1, GRATE2, GRATE3, GRATE4, GRATE5, GRATE6 DOUBLE PRECISION TLAM2, TLAM3, TLAM4, TLAM5, TLAM6 * THE RATIO OF LAMBDA TO PAVAL DOUBLE PRECISION LRATIO * USED TO STORE THE INVERSE OF THE 2X2 PIVOT DOUBLE PRECISION A11, A21, A22 * THE DETERMINANT OF THE 2X2 PIVOT DOUBLE PRECISION DETERM * MU AND LAMBDA FROM THE JONES-PATRICK PAPER DOUBLE PRECISION MU, LAMBDA(6,1) * THE PIVOTING DECISION CRITERION FROM BUNCH-KAUFMAN DOUBLE PRECISION CRIT * HOW MANY STEPS TO SKIP INTEGER NSKIP * SOME COUNTING VARIABLES TO GATHER INFORMATION (WILL GO) * A LOCAL VARIABLE FOR LOOP UNROLLING INTEGER IADDR INTEGER ME, NPROCS * TEMPORARY LEVEL OF LOOP UNROLLING INTEGER TNROLL * COMMUNICATION VARIABLE * INTEGER FUNCTIONS INTEGER IGETMX, LOCAL2 EXTERNAL IGETMX, LOCAL2, L2TMON, L2TMOF, CPYVEC * LINE NUMBER CONVENTIONS * 0-99: TEMPORARY CODE * 100-102: MAIN LOOP AND END OF MAIN LOOP * 110-199: READING IN MATRIX ELEMENTS THE 1ST TIME * 200-299: PIVOT SIZE DETERMINATION * 300-399: SEARCH FOR MU * 400-449: 1X1 PIVOT * 450-549: CHANGE AND READ MATRIX FOR A B-K 2X2 PIVOT * 550-599: B-K 2X2 PIVOT * 600-649: 2X2 PIVOT STEP * 650-699: 3X3 PIVOT STEP * 700-749: 4X4 PIVOT STEP * 1000-1049: 5X5 PIVOT STEP * 1050-1099: 6X6 PIVOT STEP * 750-799: EMPTY PIVOT STEP * 800-899: FORWARD, BACK, AND DIAGONAL SOLVES * 8000-8999: MORE BACK SOLVE * 1500-1549: REGROUPING AFTER FACTORIZATION * 4000-4600: PARALLEL SYNC. * 5000-6000: FIND NEXT LAMBDA IN PARALLEL * FACTOR THE MATRIX IF IT HASN'T BEEN DONE BEFORE IF ((INFLAG.EQ.0).OR.(INFLAG.EQ.6)) THEN ME = 1 NPROCS = 1 UPTO = 0 * TIME THE FACTORIZATION INTERNALLY IF (DEBUG.EQ.-2) THEN CALL L2TMON() ENDIF NEXT = 1 CRIT = 0.525D0 NSKIP = 0 SLVL1 = 1.0D0 + (1.0D0/CRIT) SLVL2 = SLVL1*SLVL1 SLVL3 = SLVL1*SLVL2 SLVL4 = SLVL2*SLVL2 SLVL5 = SLVL1*SLVL4 SLVL6 = SLVL3*SLVL3 GOTO (1,2,3,4,5,6) UNROLL PRINT *,'ERROR: ONLY LEVELS 1-6 ARE ALLOWED IN LOOP UNROLLING' STOP 1 SLVLM = SLVL1 GOTO 7 2 SLVLM = SLVL2 GOTO 7 3 SLVLM = SLVL3 GOTO 7 4 SLVLM = SLVL4 GOTO 7 5 SLVLM = SLVL5 GOTO 7 6 SLVLM = SLVL6 GOTO 7 7 CONTINUE INERT = 0 DO 20 I = 1, UNROLL LAMBDA(I,ME) = 0.0 20 CONTINUE TNROLL = UNROLL * START THE BUNCH FACTORIZATION LOOP DO 100 L = 1, N * SKIP IF A 2X2 PIVOT WAS DONE LAST STEP IF (NSKIP.GT.0) THEN NSKIP = NSKIP-1 GOTO 102 ENDIF * DETERMINE WHAT ROWS NEED TO BE READ IN * GET ALL THE NECESSARY ROWS FROM K-SIGMA*M RSTART = RWSTRT(L) DO 110 I = L+1, MIN0(L+UNROLL-1,N) RSTART = MAX0(RSTART,RWSTRT(I)) 110 CONTINUE IF ((RSTART-L.LT.12).OR.(UPTO-L.LT.UNROLL)) THEN RSTART = RWSTRT(L) TNROLL = 1 ELSE TNROLL = UNROLL ENDIF * SET ALL THE ROWPTR'S (WE DO THIS HERE FOR PARALLELISM) DO 115 I = UPTO+1, RSTART ROWPTR(I) = NEXT NEXT = NEXT + I-L+1 ROWPTR(I+1) = NEXT IF (NEXT.GT.MXSIZE) THEN PRINT *,'ERROR: B-K FACTORIZATION OUT OF SPACE', C ' AT STEP ',L,' OF ',N,NEXT,MXSIZE,INERT STOP ENDIF 115 CONTINUE DO 120 I = UPTO+1, RSTART DO 130 J = ROWPTR(I),ROWPTR(I+1)-1 A(J) = 0.0D0 130 CONTINUE START = ROWPTR(I)+I DO 140 J = KRP(I), KRP(I+1)-1 A(START-KCP(J)) = K(J) 140 CONTINUE CNORECUR DO 150 J = MRP(I), MRP(I+1)-1 A(START-MCP(J)) = A(START-MCP(J)) - SIGMA*M(J) 150 CONTINUE * FIND THE LAMBDA'S IN THE NEW DATA VL(I) = A(ROWPTR(I)+I-L) IF (LAMBDA(1,ME).LT.ABS(VL(I))) THEN TR(ME) = I LAMBDA(1,ME) = ABS(VL(I)) ENDIF IF (TNROLL.GT.1) THEN VL2(I) = A(ROWPTR(I)+I-(L+1)) LAMBDA(2,ME) = MAX(LAMBDA(2,ME),ABS(VL2(I))) ENDIF 120 CONTINUE UPTO = MAX0(UPTO,RSTART) * DETERMINE WHAT SIZE PIVOT TO USE * FIRST, FIND LAMBDA OF COLUMN L, AND PUT COLUMN L IN V R = TR(1) DO 200 I = 2, NPROCS IF (LAMBDA(1,I).GT.LAMBDA(1,1)) THEN R = TR(I) LAMBDA(1,1) = LAMBDA(1,I) ENDIF 200 CONTINUE IF (TNROLL.GT.1) THEN DO 210 I = 2, NPROCS LAMBDA(2,1) = MAX(LAMBDA(2,1),LAMBDA(2,I)) 210 CONTINUE ENDIF * IF LAMBDA IS NOT ZERO THEN TRY TO DETERMINE WHAT * SIZE PIVOT TO USE IF (LAMBDA(1,1).NE.0.0D0) THEN * CHECK THE GROWTH RATE FOR A 1X1 PAVAL = A(ROWPTR(L)) IF (PAVAL.EQ.0.0D0) GOTO 300 PCVAL = A(ROWPTR(L+1)) LRATIO = ABS(LAMBDA(1,1)/PAVAL) GRATE1 = 1.0D0 + LRATIO * GROWTH RATE OF 1X1 IS TOO BIG, DO A MU SEARCH IF (GRATE1.GT.SLVLM) GOTO 300 * DO A STANDARD B-K PIVOT CHECK IF PCVAL IS 0 * OR IF WE ARE NEAR THE END OF THE FACTORIZATION IF ((PCVAL.EQ.0.0D0).OR.(UPTO-L.LT.12).OR. C (TNROLL.EQ.1)) GOTO 264 * FIND THE 2X2 GROWTH RATE TBVAL = A(ROWPTR(L+1)+1) PBVAL = TBVAL/PAVAL PCVAL = PCVAL - (PBVAL*TBVAL) TLAM2 = ABS(PBVAL)*LAMBDA(1,1)+LAMBDA(2,1) GRATE2 = GRATE1*(1.0D0+TLAM2/ABS(PCVAL)) * DO A PRELIMINARY CHECK ON GRATE2, IF IT IS TOO LARGE * ALREADY FOR A 6X6 THEN JUST GO THE 1X1 TEST IF ((GRATE2.GT.SLVLM).OR.(LAMBDA(2,1).EQ.0.0D0)) C GOTO 264 GOTO 263 263 CONTINUE * TRY TO DO A 2X2 PIVOT IF (GRATE2.LE.SLVL2) THEN GOTO 600 ENDIF 264 CONTINUE * TRY TO USE A(L,L) FOR A 1 X 1 PIVOT IF (GRATE1.LE.SLVL1) THEN GOTO 400 ENDIF * END OF PIVOT SIZE DETERMINATION * OTHERWISE WE FALL THROUGH AND SEARCH FOR MU ELSE * THIS ELSE STATEMENT CORRESPONDS WITH LAMBDA=0 * 0X0 PIVOT GOTO 750 * END OF 0 PIVOT ENDIF 300 CONTINUE * MU SEARCH * DIDN'T PASS ABOVE TESTS, NOW WE MUST FIND SIGMA * CALLED MU HERE TO AVOID CONFUSION WITH OTHER SIGMA * FIRST SEARCH ROW R, THEN COLUMN R TO DETERMINE MU MU = 0.0D0 DO 310 I = ROWPTR(R), ROWPTR(R)+(R-L)-1 MU = MAX(ABS(A(I)),MU) 310 CONTINUE TMPMAX = MAX0(UPTO,RWSTRT(R)) IF (UPTO.LT.RWSTRT(R)) THEN DO 320 I = UPTO+1, TMPMAX V(I) = 0.0D0 320 CONTINUE DO 330 I = CKRP(R), CKRP(R+1)-1 V(CKCP(I)) = CKOFF(I) 330 CONTINUE CNORECUR DO 340 I = CMRP(R), CMRP(R+1)-1 J = CMCP(I) V(J) = V(J) - SIGMA*CMOFF(I) 340 CONTINUE ENDIF DO 350 I = R+1, UPTO V(I) = A(ROWPTR(I)+I-R) 350 CONTINUE DO 360 I = R+1, TMPMAX MU = MAX(MU,ABS(V(I))) 360 CONTINUE * END OF MU SEARCH * NOW CHECK TO SEE IF WE CAN STILL USE A(L,L) AS THE PIVOT IF (CRIT*LAMBDA(1,1)*LAMBDA(1,1).LE. C MU*ABS(A(ROWPTR(L)))) THEN * GO AHEAD AND USE A(L,L) AND DO THE 1X1 PIVOT GOTO 400 ELSE GOTO 450 ENDIF 400 CONTINUE * THE 1X1 PIVOT IUPTO(L) = UPTO DIAG(L) = 1.0D0/A(ROWPTR(L)) IF (DIAG(L).LT.0.0D0) THEN INERT = INERT + 1 ENDIF A(ROWPTR(L)) = 1.0D0 SUBDAG(L+1) = 0.0D0 P(L) = 3 AT = DIAG(L) * COMPUTE THE PIVOT ROW CNORECUR DO 410 I = L+1, UPTO V(I) = VL(I) 410 CONTINUE * UPDATE THE ROWS BELOW DO 420 I = L+1, UPTO T = V(I) * AT START = ROWPTR(I)+I A(START-L) = T DO 430 J = L+1, I A(START-J) = A(START-J) - T*V(J) 430 CONTINUE 420 CONTINUE NSKIP = 0 GOTO 101 * END OF 1X1 PIVOT 450 CONTINUE * CHANGE MATRIX FOR BUNCH-KAUFMAN'S 2X2 PIVOT * DIDN'T PASS ABOVE TESTS, SO NOW WE MUST DO A 2X2 * FIRST, INTERCHANGE ROWS AND COLUMNS L+1 AND R * READ THE NECESSARY ROWS FROM K-SIGMA*M * AND ALLOW FOR THE FILL-IN IF (DEBUG.GT.1) THEN PRINT *,'FILLING IN ',RWSTRT(R)-UPTO,' AT ',L,R,UPTO ENDIF * IF WE DIDN'T EVEN TRY TO LOOK AT A 2X2 THEN FILL VL2 IF (TNROLL.EQ.1) THEN DO 453 I = L+1, UPTO VL2(I) = A(ROWPTR(I)+I-(L+1)) 453 CONTINUE ENDIF * GET ALL THE NECESSARY ROWS FROM K-SIGMA*M DO 454 I = UPTO+1, RWSTRT(R) ROWPTR(I) = NEXT NEXT = NEXT + I-L+1 ROWPTR(I+1) = NEXT IF (NEXT.GT.MXSIZE) THEN PRINT *,'ERROR: B-K FACTORIZATION OUT OF SPACE', C ' AT STEP ',L,' OF ',N,NEXT,MXSIZE,INERT STOP ENDIF 454 CONTINUE DO 455 I = UPTO+1, RWSTRT(R) DO 460 J = ROWPTR(I),ROWPTR(I+1)-1 A(J) = 0.0D0 460 CONTINUE START = ROWPTR(I)+I DO 470 J = KRP(I), KRP(I+1)-1 A(START-KCP(J)) = K(J) 470 CONTINUE CNORECUR DO 480 J = MRP(I), MRP(I+1)-1 A(START-MCP(J)) = A(START-MCP(J)) - SIGMA*M(J) 480 CONTINUE * SWITCH COLUMN L+1 AND R IN THIS NEW ROW (I) * PUTTING COLUMN L+1 INTO V2 V2(I) = A(START-R) A(START-R) = A(START-(L+1)) 455 CONTINUE * TACK ON ANY NECESSARY TRAILING 0'S TO V DO 490 I = UPTO+1, RWSTRT(R) VL(I) = 0.0D0 490 CONTINUE * SWITCH COLUMN L+1 AND R IN ROWS FROM R+1 TO UPTO CNORECUR DO 500 I = R+1, UPTO IADDR = ROWPTR(I)+I V2(I) = A(IADDR-R) A(IADDR-R) = VL2(I) 500 CONTINUE UPTO = MAX0(UPTO,RWSTRT(R)) IADDR = ROWPTR(R)+R * ROW AND COLUMN SWITCH IN ROWS L+2 TO R-1 CNORECUR DO 510 I = L+2, R-1 V2(I) = A(IADDR-I) A(IADDR-I) = VL2(I) 510 CONTINUE * SWITCH A(L+1,L+1) WITH A(R,R) V2(L+1) = A(ROWPTR(R)) A(ROWPTR(R)) = VL2(L+1) * FIX UP THE V VECTOR SO THAT IT CONTAINS COLUMN L * EFFECTIVELY SWITCH A(L+1,L) WITH A(R,L) T = VL(L+1) VL(L+1) = VL(R) VL(R) = T * END OF MATRIX CHANGES * BUNCH-KAUFMAN 2X2 PIVOT STEP * LOAD V2(R) V2(R) = VL2(R) INERT = INERT + 1 P(L) = 2 P(L+1) = -R T1 = VL(L) T2 = V2(L+1) T3 = VL(L+1) IUPTO(L) = UPTO IUPTO(L+1) = UPTO * COMPUTE THE INVERSE OF THE 2X2 PIVOT DETERM = ((((T1*T2)/T3)-T3)*T3) A11 = T2/DETERM A22 = T1/DETERM A21 = -T3/DETERM DIAG(L) = A11 DIAG(L+1) = A22 SUBDAG(L+1) = A21 SUBDAG(L+2) = 0.0D0 A(ROWPTR(L)) = 1.0D0 A(ROWPTR(L+1)) = 1.0D0 A(ROWPTR(L+1)+1) = 0.0D0 * COMPUTE THE TWO PIVOT ROWS CNORECUR DO 570 I = L+2, UPTO V(I) = VL(I) 570 CONTINUE * UPDATE THE ROWS BELOW CNORECUR DO 579 I = L+2, UPTO START = ROWPTR(I)+I VL(I) = V(I)*A11 + V2(I)*A21 A(START-L) = VL(I) VL2(I) = V(I)*A21 + V2(I)*A22 A(START-(L+1)) = VL2(I) 579 CONTINUE DO 580 I = L+2, UPTO START = ROWPTR(I)+I T = VL(I) T1 = VL2(I) CNORECUR DO 590 J = L+2, I A(START-J) = A(START-J) - T*V(J) - T1*V2(J) 590 CONTINUE 580 CONTINUE NSKIP = 1 GOTO 101 * END OF BUNCH-KAUFMAN 2X2 PIVOT STEP 600 CONTINUE * 2X2 PIVOT IUPTO(L) = UPTO IUPTO(L+1) = UPTO P(L) = 4 P(L+1) = 4 DIAG(L) = 1.0D0/PAVAL DIAG(L+1) = 1.0D0/PCVAL IF (PAVAL.LT.0.0D0) THEN INERT = INERT + 1 ENDIF IF (PCVAL.LT.0.0D0) THEN INERT = INERT + 1 ENDIF A(ROWPTR(L)) = 1.0D0 A(ROWPTR(L+1)) = 1.0D0 SUBDAG(L+1) = 0.0D0 SUBDAG(L+2) = 0.0D0 AT = DIAG(L) AT1 = DIAG(L+1) A(ROWPTR(L+1)+1) = PBVAL CNORECUR DO 610 I = L+2, UPTO V(I) = VL(I) V2(I) = VL2(I) - VL(I)*PBVAL 610 CONTINUE DO 620 I = L+2, UPTO START = ROWPTR(I)+I T = V(I)*AT A(START-L) = T T1 = V2(I)*AT1 A(START-(L+1)) = T1 DO 630 J = L+2, I A(START-J) = A(START-J) - T*V(J) - T1*V2(J) 630 CONTINUE 620 CONTINUE NSKIP = 1 GOTO 101 * END OF 2X2 PIVOT STEP 650 CONTINUE * END OF 3X3 PIVOT STEP 700 CONTINUE * END OF 4X4 PIVOT STEP 1000 CONTINUE * END OF 5X5 PIVOT STEP 1050 CONTINUE * END OF 6X6 PIVOT STEP 750 CONTINUE * EMPTY PIVOT STEP DIAG(L) = 1.0D0 / A(ROWPTR(L)) IF (DIAG(L).LT.0.0D0) THEN INERT = INERT + 1 ENDIF SUBDAG(L+1) = 0.0D0 A(ROWPTR(L)) = 1.0D0 P(L) = 1 IUPTO(L) = L GOTO 101 * END OF EMPTY PIVOT STEP 101 CONTINUE START = L + NSKIP + 1 DO 5000 I = 1, UNROLL LAMBDA(I,ME) = 0.0 5000 CONTINUE TNROLL = MAX0(0,UPTO-START) TNROLL = MIN0(UNROLL,TNROLL) IF (TNROLL.EQ.0) GOTO 102 * FIND NEXT LAMBDA(1) NSTART = START + (ME-1) IF (ME.EQ.1) THEN VL(NSTART) = A(ROWPTR(NSTART)) NSTART = NSTART + NPROCS ENDIF DO 5010 I = NSTART, UPTO, NPROCS VL(I) = A(ROWPTR(I)+I-START) 5010 CONTINUE TR(ME) = IGETMX(NSTART,UPTO,NPROCS,VL,LAMBDA(1,ME)) GOTO (102,5100,5200,5300,5400,5500) TNROLL 5100 CONTINUE IF (NSTART.EQ.START+1) THEN VL2(NSTART) = A(ROWPTR(NSTART)) NSTART = NSTART + NPROCS ENDIF DO 5110 I = NSTART, UPTO, NPROCS VL2(I) = A(ROWPTR(I)+I-(START+1)) LAMBDA(2,ME) = MAX(LAMBDA(2,ME),ABS(VL2(I))) 5110 CONTINUE GOTO 102 5200 CONTINUE 5300 CONTINUE 5400 CONTINUE 5500 CONTINUE 102 CONTINUE 100 CONTINUE * TIME FACTORIZATION INTERNALLY IF (DEBUG.EQ.-2) THEN CALL L2TMOF() PRINT *,'TIME FOR FACTORIZATION = ',LOCAL2() ENDIF IF (DEBUG.GT.1) THEN PRINT *,'USED ',NEXT-1,' OUT OF ',MXSIZE PRINT *,'INERTIA = ',INERT,' WITH SIGMA = ',SIGMA ENDIF * COMBINE EQUAL LENGTH ROWS FOR THE SOLVES I = 1 ICMAX = 10 * GOTO 1598 1500 CONTINUE GOTO (1501,1502,1503,1504,1505,1506,1507,1508,1529,1530) P(I) * EMPTY DON'T COMBINE 1501 CONTINUE I = I + 1 GOTO 1509 * B-K DON'T COMBINE 1502 CONTINUE I = I + 2 GOTO 1509 * 1X1 1503 CONTINUE IF ((IUPTO(I).EQ.IUPTO(I+1)).AND.(P(I+1).LT.ICMAX).AND. C (P(I+1).GE.3)) THEN P(I) = P(I) + (P(I+1)-2) CDIR$ NOVECTOR DO 1510 J = I+1, I+P(I)-3 P(J) = P(I) 1510 CONTINUE CDIR$ VECTOR ELSE I = I + 1 ENDIF GOTO 1509 * 2X2 1504 CONTINUE IF ((IUPTO(I).EQ.IUPTO(I+2)).AND.(P(I+2).LT.ICMAX-1).AND. C (P(I+2).GE.3)) THEN P(I) = P(I) + (P(I+2)-2) CDIR$ NOVECTOR DO 1511 J = I+1, I+P(I)-3 P(J) = P(I) 1511 CONTINUE CDIR$ VECTOR ELSE I = I + 2 ENDIF GOTO 1509 * 3X3 1505 CONTINUE IF ((IUPTO(I).EQ.IUPTO(I+3)).AND.(P(I+3).LT.ICMAX-2).AND. C (P(I+3).GE.3)) THEN P(I) = P(I) + (P(I+3)-2) CDIR$ NOVECTOR DO 1512 J = I+1, I+P(I)-3 P(J) = P(I) 1512 CONTINUE CDIR$ VECTOR ELSE I = I + 3 ENDIF GOTO 1509 * 4X4 1506 CONTINUE IF ((IUPTO(I).EQ.IUPTO(I+4)).AND.(P(I+4).LT.ICMAX-3).AND. C (P(I+4).GE.3)) THEN P(I) = P(I) + (P(I+4)-2) CDIR$ NOVECTOR DO 1513 J = I+1, I+P(I)-3 P(J) = P(I) 1513 CONTINUE CDIR$ VECTOR ELSE I = I + 4 ENDIF GOTO 1509 * 5X5 1507 CONTINUE IF ((IUPTO(I).EQ.IUPTO(I+5)).AND.(P(I+5).LT.ICMAX-4).AND. C (P(I+5).GE.3)) THEN P(I) = P(I) + (P(I+5)-2) CDIR$ NOVECTOR DO 1514 J = I+1, I+P(I)-3 P(J) = P(I) 1514 CONTINUE CDIR$ VECTOR ELSE I = I + 5 ENDIF GOTO 1509 * 6X6 1508 CONTINUE IF ((IUPTO(I).EQ.IUPTO(I+6)).AND.(P(I+6).LT.ICMAX-5).AND. C (P(I+6).GE.3)) THEN P(I) = P(I) + (P(I+6)-2) CDIR$ NOVECTOR DO 1515 J = I+1, I+P(I)-3 P(J) = P(I) 1515 CONTINUE CDIR$ VECTOR ELSE I = I + 6 ENDIF GOTO 1509 * 7X7 1529 CONTINUE IF ((IUPTO(I).EQ.IUPTO(I+7)).AND.(P(I+7).LT.ICMAX-6).AND. C (P(I+7).GE.3)) THEN P(I) = P(I) + (P(I+7)-2) CDIR$ NOVECTOR DO 1516 J = I+1, I+P(I)-3 P(J) = P(I) 1516 CONTINUE CDIR$ VECTOR ELSE I = I + 7 ENDIF GOTO 1509 * 8X8 1530 CONTINUE IF ((IUPTO(I).EQ.IUPTO(I+8)).AND.(P(I+8).LT.ICMAX-7).AND. C (P(I+8).GE.3)) THEN P(I) = P(I) + (P(I+8)-2) CDIR$ NOVECTOR DO 1517 J = I+1, I+P(I)-3 P(J) = P(I) 1517 CONTINUE CDIR$ VECTOR ELSE I = I + 8 ENDIF GOTO 1509 1509 CONTINUE IF (I.LT.N-(ICMAX-2)) GOTO 1500 *1598 CONTINUE IF (DEBUG.GE.2) THEN IT1 = 0 IT2 = 0 IT3 = 0 IT4 = 0 IT5 = 0 IT6 = 0 IT7 = 0 IT8 = 0 DO 1599 I = 1, N IF (P(I).EQ.3) THEN IT1 = IT1 + 1 ELSE IF (P(I).EQ.4) THEN IT2 = IT2 + 1 ELSE IF (P(I).EQ.5) THEN IT3 = IT3 + 1 ELSE IF (P(I).EQ.6) THEN IT4 = IT4 + 1 ELSE IF (P(I).EQ.7) THEN IT5 = IT5 + 1 ELSE IF (P(I).EQ.8) THEN IT6 = IT6 + 1 ELSE IF (P(I).EQ.9) THEN IT7 = IT7 + 1 ELSE IF (P(I).EQ.10) THEN IT8 = IT8 + 1 ENDIF 1599 CONTINUE PRINT *,'COMB 1 = ',IT1 PRINT *,'COMB 2 = ',IT2/2 PRINT *,'COMB 3 = ',IT3/3 PRINT *,'COMB 4 = ',IT4/4 PRINT *,'COMB 5 = ',IT5/5 PRINT *,'COMB 6 = ',IT6/6 PRINT *,'COMB 7 = ',IT7/7 PRINT *,'COMB 8 = ',IT8/8 ENDIF ENDIF IF (INFLAG.EQ.6) THEN RETURN ENDIF * TIME THE FORWARD SOLVE INTERNALLY IF (DEBUG.EQ.-2) THEN CALL L2TMON() ENDIF CALL CPYVEC(N,V,B) * SOLVE P(T)L V = B * DO 810 I = 1, N I = 1 810 CONTINUE IF (I.GT.N) GOTO 849 GOTO (811,812,813,814,815,816,817,818,819,820) P(I) 811 CONTINUE * A 0X0 PIVOT GOTO 831 812 CONTINUE * A B-K 2X2 PIVOT L = -P(I+1) T = V(L) V(L) = V(I+1) V(I+1) = T - A(ROWPTR(I+1)+1)*V(I) GOTO 834 813 CONTINUE * A 1X1 PIVOT GOTO 833 814 CONTINUE * A 2X2 PIVOT V(I+1) = V(I+1) - A(ROWPTR(I+1)+1)*V(I) GOTO 834 815 CONTINUE * A 3X3 PIVOT V(I+1) = V(I+1) - A(ROWPTR(I+1)+1)*V(I) V(I+2) = V(I+2) - A(ROWPTR(I+2)+2)*V(I) - C A(ROWPTR(I+2)+1)*V(I+1) GOTO 835 816 CONTINUE * A 4X4 PIVOT V(I+1) = V(I+1) - A(ROWPTR(I+1)+1)*V(I) V(I+2) = V(I+2) - A(ROWPTR(I+2)+2)*V(I) - C A(ROWPTR(I+2)+1)*V(I+1) V(I+3) = V(I+3) - A(ROWPTR(I+3)+3)*V(I) - C A(ROWPTR(I+3)+2)*V(I+1) - A(ROWPTR(I+3)+1)*V(I+2) GOTO 836 817 CONTINUE * A 5X5 PIVOT V(I+1) = V(I+1) - A(ROWPTR(I+1)+1)*V(I) V(I+2) = V(I+2) - A(ROWPTR(I+2)+2)*V(I) - C A(ROWPTR(I+2)+1)*V(I+1) V(I+3) = V(I+3) - A(ROWPTR(I+3)+3)*V(I) - C A(ROWPTR(I+3)+2)*V(I+1) - A(ROWPTR(I+3)+1)*V(I+2) V(I+4) = V(I+4) - A(ROWPTR(I+4)+4)*V(I) - C A(ROWPTR(I+4)+3)*V(I+1) - A(ROWPTR(I+4)+2)*V(I+2) - C A(ROWPTR(I+4)+1)*V(I+3) GOTO 837 818 CONTINUE * A 6X6 PIVOT V(I+1) = V(I+1) - A(ROWPTR(I+1)+1)*V(I) V(I+2) = V(I+2) - A(ROWPTR(I+2)+2)*V(I) - C A(ROWPTR(I+2)+1)*V(I+1) V(I+3) = V(I+3) - A(ROWPTR(I+3)+3)*V(I) - C A(ROWPTR(I+3)+2)*V(I+1) - A(ROWPTR(I+3)+1)*V(I+2) V(I+4) = V(I+4) - A(ROWPTR(I+4)+4)*V(I) - C A(ROWPTR(I+4)+3)*V(I+1) - A(ROWPTR(I+4)+2)*V(I+2) - C A(ROWPTR(I+4)+1)*V(I+3) V(I+5) = V(I+5) - A(ROWPTR(I+5)+5)*V(I) - C A(ROWPTR(I+5)+4)*V(I+1) - A(ROWPTR(I+5)+3)*V(I+2) - C A(ROWPTR(I+5)+2)*V(I+3) - A(ROWPTR(I+5)+1)*V(I+4) GOTO 838 819 CONTINUE * A 7X7 PIVOT V(I+1) = V(I+1) - A(ROWPTR(I+1)+1)*V(I) V(I+2) = V(I+2) - A(ROWPTR(I+2)+2)*V(I) - C A(ROWPTR(I+2)+1)*V(I+1) V(I+3) = V(I+3) - A(ROWPTR(I+3)+3)*V(I) - C A(ROWPTR(I+3)+2)*V(I+1) - A(ROWPTR(I+3)+1)*V(I+2) V(I+4) = V(I+4) - A(ROWPTR(I+4)+4)*V(I) - C A(ROWPTR(I+4)+3)*V(I+1) - A(ROWPTR(I+4)+2)*V(I+2) - C A(ROWPTR(I+4)+1)*V(I+3) V(I+5) = V(I+5) - A(ROWPTR(I+5)+5)*V(I) - C A(ROWPTR(I+5)+4)*V(I+1) - A(ROWPTR(I+5)+3)*V(I+2) - C A(ROWPTR(I+5)+2)*V(I+3) - A(ROWPTR(I+5)+1)*V(I+4) V(I+6) = V(I+6) - A(ROWPTR(I+6)+6)*V(I) - C A(ROWPTR(I+6)+5)*V(I+1) - A(ROWPTR(I+6)+4)*V(I+2) - C A(ROWPTR(I+6)+3)*V(I+3) - A(ROWPTR(I+6)+2)*V(I+4) - C A(ROWPTR(I+6)+1)*V(I+5) GOTO 839 820 CONTINUE * A 8X8 PIVOT V(I+1) = V(I+1) - A(ROWPTR(I+1)+1)*V(I) V(I+2) = V(I+2) - A(ROWPTR(I+2)+2)*V(I) - C A(ROWPTR(I+2)+1)*V(I+1) V(I+3) = V(I+3) - A(ROWPTR(I+3)+3)*V(I) - C A(ROWPTR(I+3)+2)*V(I+1) - A(ROWPTR(I+3)+1)*V(I+2) V(I+4) = V(I+4) - A(ROWPTR(I+4)+4)*V(I) - C A(ROWPTR(I+4)+3)*V(I+1) - A(ROWPTR(I+4)+2)*V(I+2) - C A(ROWPTR(I+4)+1)*V(I+3) V(I+5) = V(I+5) - A(ROWPTR(I+5)+5)*V(I) - C A(ROWPTR(I+5)+4)*V(I+1) - A(ROWPTR(I+5)+3)*V(I+2) - C A(ROWPTR(I+5)+2)*V(I+3) - A(ROWPTR(I+5)+1)*V(I+4) V(I+6) = V(I+6) - A(ROWPTR(I+6)+6)*V(I) - C A(ROWPTR(I+6)+5)*V(I+1) - A(ROWPTR(I+6)+4)*V(I+2) - C A(ROWPTR(I+6)+3)*V(I+3) - A(ROWPTR(I+6)+2)*V(I+4) - C A(ROWPTR(I+6)+1)*V(I+5) V(I+7) = V(I+7) - A(ROWPTR(I+7)+7)*V(I) - C A(ROWPTR(I+7)+6)*V(I+1) - A(ROWPTR(I+7)+5)*V(I+2) - C A(ROWPTR(I+7)+4)*V(I+3) - A(ROWPTR(I+7)+3)*V(I+4) - C A(ROWPTR(I+7)+2)*V(I+5) - A(ROWPTR(I+7)+1)*V(I+6) GOTO 840 831 CONTINUE * A 0X0 PIVOT I = I + 1 GOTO 810 833 CONTINUE * A 1X1 PIVOT CNORECUR DO 841 J = I+1, IUPTO(I) V(J) = V(J) - A(ROWPTR(J)+J-I)*V(I) 841 CONTINUE I = I + 1 GOTO 810 834 CONTINUE * A 2X2 PIVOT CNORECUR DO 842 J = I+2, IUPTO(I) IADDR = ROWPTR(J)+J-I V(J) = V(J) - A(IADDR)*V(I) - A(IADDR-1)*V(I+1) 842 CONTINUE I = I + 2 GOTO 810 835 CONTINUE * A 3X3 PIVOT CNORECUR DO 843 J = I+3, IUPTO(I) IADDR = ROWPTR(J)+J-I V(J) = V(J) - A(IADDR)*V(I) - A(IADDR-1)*V(I+1) - C A(IADDR-2)*V(I+2) 843 CONTINUE I = I + 3 GOTO 810 836 CONTINUE * A 4X4 PIVOT CNORECUR DO 844 J = I+4, IUPTO(I) IADDR = ROWPTR(J)+J-I V(J) = V(J) - A(IADDR)*V(I) - A(IADDR-1)*V(I+1) - C A(IADDR-2)*V(I+2) - A(IADDR-3)*V(I+3) 844 CONTINUE I = I + 4 GOTO 810 837 CONTINUE * A 5X5 PIVOT CNORECUR DO 845 J = I+5, IUPTO(I) IADDR = ROWPTR(J)+J-I V(J) = V(J) - A(IADDR)*V(I) - A(IADDR-1)*V(I+1) - C A(IADDR-2)*V(I+2) - A(IADDR-3)*V(I+3) - A(IADDR-4)*V(I+4) 845 CONTINUE I = I + 5 GOTO 810 838 CONTINUE * A 6X6 PIVOT CNORECUR DO 846 J = I+6, IUPTO(I) IADDR = ROWPTR(J)+J-I V(J) = V(J) - A(IADDR)*V(I) - A(IADDR-1)*V(I+1) - C A(IADDR-2)*V(I+2) - C A(IADDR-3)*V(I+3) - A(IADDR-4)*V(I+4) - A(IADDR-5)*V(I+5) 846 CONTINUE I = I + 6 GOTO 810 839 CONTINUE * A 7X7 PIVOT CNORECUR DO 847 J = I+7, IUPTO(I) IADDR = ROWPTR(J)+J-I V(J) = V(J) - A(IADDR)*V(I) - A(IADDR-1)*V(I+1) - C A(IADDR-2)*V(I+2) - C A(IADDR-3)*V(I+3) - A(IADDR-4)*V(I+4) - A(IADDR-5)*V(I+5) - C A(IADDR-6)*V(I+6) 847 CONTINUE I = I + 7 GOTO 810 840 CONTINUE * A 8X8 PIVOT CNORECUR DO 848 J = I+8, IUPTO(I) IADDR = ROWPTR(J)+J-I V(J) = V(J) - A(IADDR)*V(I) - A(IADDR-1)*V(I+1) - C A(IADDR-2)*V(I+2) - C A(IADDR-3)*V(I+3) - A(IADDR-4)*V(I+4) - A(IADDR-5)*V(I+5) - C A(IADDR-6)*V(I+6) - A(IADDR-7)*V(I+7) 848 CONTINUE I = I + 8 GOTO 810 * END OF FORWARD SOLVE 849 CONTINUE * TIME FORWARD SOLVE INTERNALLY IF (DEBUG.EQ.-2) THEN CALL L2TMOF() PRINT *,'TIME FOR FORWARD SOLVE = ',LOCAL2() ENDIF * SOLVE D U = V X(1) = DIAG(1)*V(1)+SUBDAG(2)*V(2) X(N) = DIAG(N)*V(N)+SUBDAG(N)*V(N-1) DO 850 I = 2, N-1 X(I) = DIAG(I)*V(I)+SUBDAG(I)*V(I-1)+SUBDAG(I+1)*V(I+1) 850 CONTINUE * TIME THE BACKWARD SOLVE INTERNALLY IF (DEBUG.EQ.-2) THEN CALL L2TMON() ENDIF * SOLVE L(T)*P X = V * DO 860 I = N, 1, -1 I = N 860 CONTINUE IF (I.LE.0) GOTO 880 IF (P(I).LE.1) THEN IF (P(I).LT.0) THEN NUMROW = 2 ELSE I = I - 1 GOTO 860 ENDIF ELSE NUMROW = P(I) - 2 ENDIF IF (NUMROW.GE.NPROCS) THEN MYNUM = NUMROW / NPROCS ISTART = MYNUM * (ME-1) IF (ME.LE.MOD(NUMROW,NPROCS)) THEN MYNUM = MYNUM + 1 ISTART = ISTART + (ME-1) ELSE ISTART = ISTART + MOD(NUMROW,NPROCS) ENDIF GOTO (861,862,863,864,865,866,867,868) MYNUM GOTO 879 ELSE ISTART = MOD((ME-1),NUMROW) IMULT = NPROCS / NUMROW IF (ISTART.LT.MOD(NPROCS,NUMROW)) THEN IMULT = IMULT + 1 ENDIF IOFF = (ME-1) / NUMROW IF (IMULT.EQ.1) THEN GOTO 861 ELSE GOTO 869 ENDIF ENDIF 861 CONTINUE * A 1X1 PIVOT T = X(I-ISTART) DO 871 J = I+1, IUPTO(I) T = T - A(ROWPTR(J)+J-I+ISTART)*X(J) 871 CONTINUE X(I-ISTART) = T GOTO 879 * END OF 1X1 862 CONTINUE * A NORMAL 2X2 PIVOT T = X(I-ISTART) T1 = X(I-ISTART-1) DO 872 J = I+1, IUPTO(I) IADDR = ROWPTR(J)+J-I+ISTART TX = X(J) T = T - A(IADDR)*TX T1 = T1 - A(IADDR+1)*TX 872 CONTINUE X(I-ISTART) = T X(I-ISTART-1) = T1 GOTO 879 * END OF NORMAL 2X2 863 CONTINUE * A NORMAL 3X3 PIVOT T = X(I-ISTART) T1 = X(I-ISTART-1) T2 = X(I-ISTART-2) DO 873 J = I+1, IUPTO(I) IADDR = ROWPTR(J)+J-I+ISTART TX = X(J) T = T - A(IADDR)*TX T1 = T1 - A(IADDR+1)*TX T2 = T2 - A(IADDR+2)*TX 873 CONTINUE X(I-ISTART) = T X(I-ISTART-1) = T1 X(I-ISTART-2) = T2 GOTO 879 * END OF NORMAL 3X3 864 CONTINUE * A NORMAL 4X4 PIVOT T = X(I-ISTART) T1 = X(I-ISTART-1) T2 = X(I-ISTART-2) T3 = X(I-ISTART-3) DO 874 J = I+1, IUPTO(I) IADDR = ROWPTR(J)+J-I+ISTART TX = X(J) T = T - A(IADDR)*TX T1 = T1 - A(IADDR+1)*TX T2 = T2 - A(IADDR+2)*TX T3 = T3 - A(IADDR+3)*TX 874 CONTINUE X(I-ISTART) = T X(I-ISTART-1) = T1 X(I-ISTART-2) = T2 X(I-ISTART-3) = T3 GOTO 879 * END OF NORMAL 4X4 865 CONTINUE * A NORMAL 5X5 PIVOT T = X(I-ISTART) T1 = X(I-ISTART-1) T2 = X(I-ISTART-2) T3 = X(I-ISTART-3) T4 = X(I-ISTART-4) DO 875 J = I+1, IUPTO(I) IADDR = ROWPTR(J)+J-I+ISTART TX = X(J) T = T - A(IADDR)*TX T1 = T1 - A(IADDR+1)*TX T2 = T2 - A(IADDR+2)*TX T3 = T3 - A(IADDR+3)*TX T4 = T4 - A(IADDR+4)*TX 875 CONTINUE X(I-ISTART) = T X(I-ISTART-1) = T1 X(I-ISTART-2) = T2 X(I-ISTART-3) = T3 X(I-ISTART-4) = T4 GOTO 879 * END OF NORMAL 5X5 866 CONTINUE * A NORMAL 6X6 PIVOT T = X(I-ISTART) T1 = X(I-ISTART-1) T2 = X(I-ISTART-2) T3 = X(I-ISTART-3) T4 = X(I-ISTART-4) T5 = X(I-ISTART-5) DO 876 J = I+1, IUPTO(I) IADDR = ROWPTR(J)+J-I+ISTART TX = X(J) T = T - A(IADDR)*TX T1 = T1 - A(IADDR+1)*TX T2 = T2 - A(IADDR+2)*TX T3 = T3 - A(IADDR+3)*TX T4 = T4 - A(IADDR+4)*TX T5 = T5 - A(IADDR+5)*TX 876 CONTINUE X(I-ISTART) = T X(I-ISTART-1) = T1 X(I-ISTART-2) = T2 X(I-ISTART-3) = T3 X(I-ISTART-4) = T4 X(I-ISTART-5) = T5 GOTO 879 * END OF NORMAL 6X6 867 CONTINUE * A NORMAL 7X7 PIVOT T = X(I-ISTART) T1 = X(I-ISTART-1) T2 = X(I-ISTART-2) T3 = X(I-ISTART-3) T4 = X(I-ISTART-4) T5 = X(I-ISTART-5) T6 = X(I-ISTART-6) DO 877 J = I+1, IUPTO(I) IADDR = ROWPTR(J)+J-I+ISTART TX = X(J) T = T - A(IADDR)*TX T1 = T1 - A(IADDR+1)*TX T2 = T2 - A(IADDR+2)*TX T3 = T3 - A(IADDR+3)*TX T4 = T4 - A(IADDR+4)*TX T5 = T5 - A(IADDR+5)*TX T6 = T6 - A(IADDR+6)*TX 877 CONTINUE X(I-ISTART) = T X(I-ISTART-1) = T1 X(I-ISTART-2) = T2 X(I-ISTART-3) = T3 X(I-ISTART-4) = T4 X(I-ISTART-5) = T5 X(I-ISTART-6) = T6 GOTO 879 * END OF NORMAL 7X7 868 CONTINUE * A NORMAL 8X8 PIVOT T = X(I-ISTART) T1 = X(I-ISTART-1) T2 = X(I-ISTART-2) T3 = X(I-ISTART-3) T4 = X(I-ISTART-4) T5 = X(I-ISTART-5) T6 = X(I-ISTART-6) T7 = X(I-ISTART-7) DO 878 J = I+1, IUPTO(I) IADDR = ROWPTR(J)+J-I+ISTART TX = X(J) T = T - A(IADDR)*TX T1 = T1 - A(IADDR+1)*TX T2 = T2 - A(IADDR+2)*TX T3 = T3 - A(IADDR+3)*TX T4 = T4 - A(IADDR+4)*TX T5 = T5 - A(IADDR+5)*TX T6 = T6 - A(IADDR+6)*TX T7 = T7 - A(IADDR+7)*TX 878 CONTINUE X(I-ISTART) = T X(I-ISTART-1) = T1 X(I-ISTART-2) = T2 X(I-ISTART-3) = T3 X(I-ISTART-4) = T4 X(I-ISTART-5) = T5 X(I-ISTART-6) = T6 X(I-ISTART-7) = T7 GOTO 879 * END OF NORMAL 8X8 869 CONTINUE * A 1X1 PIVOT WITH MULTIPLE PROCESSORS T = 0.0D0 DO 889 J = I+1+IOFF, IUPTO(I), IMULT T = T + A(ROWPTR(J)+J-I+ISTART)*X(J) 889 CONTINUE X(I-ISTART) = X(I-ISTART) - T GOTO 879 * END OF 1X1 WITH MULTIPLE PROCESSORS 879 CONTINUE I = I - NUMROW II = I + NUMROW GOTO (8499,8502,8503,8504,8505,8506,8507,8508) NUMROW * 2X2 PIVOT 8502 CONTINUE X(II-1) = X(II-1) - A(ROWPTR(II)+1)*X(II) * MAKE A CORRECTION FOR B-K IF NECESSARY IF (P(II).NE.4) THEN L = -P(II) T = X(II) X(II) = X(L) X(L) = T ENDIF GOTO 8499 * 3X3 PIVOT 8503 CONTINUE X(II-1) = X(II-1) - A(ROWPTR(II)+1)*X(II) X(II-2) = X(II-2) - A(ROWPTR(II)+2)*X(II) - C A(ROWPTR(II-1)+1)*X(II-1) GOTO 8499 * 4X4 PIVOT 8504 CONTINUE X(II-1) = X(II-1) - A(ROWPTR(II)+1)*X(II) X(II-2) = X(II-2) - A(ROWPTR(II)+2)*X(II) - C A(ROWPTR(II-1)+1)*X(II-1) X(II-3) = X(II-3) - A(ROWPTR(II)+3)*X(II) - C A(ROWPTR(II-1)+2)*X(II-1) - A(ROWPTR(II-2)+1)*X(II-2) GOTO 8499 * 5X5 PIVOT 8505 CONTINUE X(II-1) = X(II-1) - A(ROWPTR(II)+1)*X(II) X(II-2) = X(II-2) - A(ROWPTR(II)+2)*X(II) - C A(ROWPTR(II-1)+1)*X(II-1) X(II-3) = X(II-3) - A(ROWPTR(II)+3)*X(II) - C A(ROWPTR(II-1)+2)*X(II-1) - A(ROWPTR(II-2)+1)*X(II-2) X(II-4) = X(II-4) - A(ROWPTR(II)+4)*X(II) - C A(ROWPTR(II-1)+3)*X(II-1) - A(ROWPTR(II-2)+2)*X(II-2) - C A(ROWPTR(II-3)+1)*X(II-3) GOTO 8499 8506 CONTINUE X(II-1) = X(II-1) - A(ROWPTR(II)+1)*X(II) X(II-2) = X(II-2) - A(ROWPTR(II)+2)*X(II) - C A(ROWPTR(II-1)+1)*X(II-1) X(II-3) = X(II-3) - A(ROWPTR(II)+3)*X(II) - C A(ROWPTR(II-1)+2)*X(II-1) - A(ROWPTR(II-2)+1)*X(II-2) X(II-4) = X(II-4) - A(ROWPTR(II)+4)*X(II) - C A(ROWPTR(II-1)+3)*X(II-1) - A(ROWPTR(II-2)+2)*X(II-2) - C A(ROWPTR(II-3)+1)*X(II-3) X(II-5) = X(II-5) - A(ROWPTR(II)+5)*X(II) - C A(ROWPTR(II-1)+4)*X(II-1) - A(ROWPTR(II-2)+3)*X(II-2) - C A(ROWPTR(II-3)+2)*X(II-3) - A(ROWPTR(II-4)+1)*X(II-4) GOTO 8499 8507 CONTINUE X(II-1) = X(II-1) - A(ROWPTR(II)+1)*X(II) X(II-2) = X(II-2) - A(ROWPTR(II)+2)*X(II) - C A(ROWPTR(II-1)+1)*X(II-1) X(II-3) = X(II-3) - A(ROWPTR(II)+3)*X(II) - C A(ROWPTR(II-1)+2)*X(II-1) - A(ROWPTR(II-2)+1)*X(II-2) X(II-4) = X(II-4) - A(ROWPTR(II)+4)*X(II) - C A(ROWPTR(II-1)+3)*X(II-1) - A(ROWPTR(II-2)+2)*X(II-2) - C A(ROWPTR(II-3)+1)*X(II-3) X(II-5) = X(II-5) - A(ROWPTR(II)+5)*X(II) - C A(ROWPTR(II-1)+4)*X(II-1) - A(ROWPTR(II-2)+3)*X(II-2) - C A(ROWPTR(II-3)+2)*X(II-3) - A(ROWPTR(II-4)+1)*X(II-4) X(II-6) = X(II-6) - A(ROWPTR(II)+6)*X(II) - C A(ROWPTR(II-1)+5)*X(II-1) - A(ROWPTR(II-2)+4)*X(II-2) - C A(ROWPTR(II-3)+3)*X(II-3) - A(ROWPTR(II-4)+2)*X(II-4) - C A(ROWPTR(II-5)+1)*X(II-5) GOTO 8499 8508 CONTINUE X(II-1) = X(II-1) - A(ROWPTR(II)+1)*X(II) X(II-2) = X(II-2) - A(ROWPTR(II)+2)*X(II) - C A(ROWPTR(II-1)+1)*X(II-1) X(II-3) = X(II-3) - A(ROWPTR(II)+3)*X(II) - C A(ROWPTR(II-1)+2)*X(II-1) - A(ROWPTR(II-2)+1)*X(II-2) X(II-4) = X(II-4) - A(ROWPTR(II)+4)*X(II) - C A(ROWPTR(II-1)+3)*X(II-1) - A(ROWPTR(II-2)+2)*X(II-2) - C A(ROWPTR(II-3)+1)*X(II-3) X(II-5) = X(II-5) - A(ROWPTR(II)+5)*X(II) - C A(ROWPTR(II-1)+4)*X(II-1) - A(ROWPTR(II-2)+3)*X(II-2) - C A(ROWPTR(II-3)+2)*X(II-3) - A(ROWPTR(II-4)+1)*X(II-4) X(II-6) = X(II-6) - A(ROWPTR(II)+6)*X(II) - C A(ROWPTR(II-1)+5)*X(II-1) - A(ROWPTR(II-2)+4)*X(II-2) - C A(ROWPTR(II-3)+3)*X(II-3) - A(ROWPTR(II-4)+2)*X(II-4) - C A(ROWPTR(II-5)+1)*X(II-5) X(II-7) = X(II-7) - A(ROWPTR(II)+7)*X(II) - C A(ROWPTR(II-1)+6)*X(II-1) - A(ROWPTR(II-2)+5)*X(II-2) - C A(ROWPTR(II-3)+4)*X(II-3) - A(ROWPTR(II-4)+3)*X(II-4) - C A(ROWPTR(II-5)+2)*X(II-5) - A(ROWPTR(II-6)+1)*X(II-6) GOTO 8499 8499 CONTINUE GOTO 860 * END OF BACKWARD SOLVE 880 CONTINUE * TIME BACKWARD SOLVE INTERNALLY IF (DEBUG.EQ.-2) THEN CALL L2TMOF() PRINT *,'TIME FOR BACKWARD SOLVE = ',LOCAL2() ENDIF RETURN END C$FORTRAN CHKFIN *********************************************************************** * LANZ SOFTWARE PACKAGE * * FILENAME: CHKFIN.MSC * * AUTHOR: MARK JONES * * PURPOSE: DETERMINES IF ANY DESIRED EIGENVALUES HAVE BEEN * * MISSED BY LOOKING AT INERTIA COUNTS, UNCONVERGED * * AND CONVERGED EIGENVALUES * *********************************************************************** SUBROUTINE CHKFIN(N,OEIGNM,OTHETA,OLDBJ,R1LEN,R1THET,R1BJ, C R2LEN,R2THET,R2BJ,ONMSGT,OSIGMA,LFTVAL,RGTVAL, C NWLEFT,NWRGHT,DEBUG,KNORM,MNORM,EPS, C ATOL,ERRCHK,NMSGT,WORKV1,SRTTHT,THTIDX, C INTLST,SIGLST,INTIDX, C FALEAD,FA,FARP,FARL,FIWRKV,FRWRKV,PROB, C NROLL,KAUX,KRP,KCP,MAUX,MRP,MCP,FASIZE,FAADDR,TMEM, C KFORM,MFORM,FACTYP,WRKPTR,RETVAL,TSTATE,CBOUND,BLEFT,BRIGHT, C BFACT,MAXDLY,INSIZE,INPTR,INADDR) * THE ORDER OF THE MATRIX INTEGER N * THE NUMBER OF EIGENVALUES IN OTHETA INTEGER OEIGNM * THE CONVERGED EIGENVALUES DOUBLE PRECISION OTHETA(*) * THE ERROR BOUNDS OF THE CONVERGED EIGENVALUES DOUBLE PRECISION OLDBJ(*) * THE NUMBER OF EIGENVALUES IN THE R1THET INTEGER R1LEN * THE GROUP OF EIGENVALUES FROM THE R1 RUN DOUBLE PRECISION R1THET(*) * THE ERROR BOUNDS FOR THE R1 RUN DOUBLE PRECISION R1BJ(*) * THE NUMBER OF EIGENVALUES IN THE R2THET INTEGER R2LEN * THE GROUP OF EIGENVALUES FROM THE R2 RUN DOUBLE PRECISION R2THET(*) * THE ERROR BOUNDS FOR THE R2 RUN DOUBLE PRECISION R2BJ(*) * THE ORIGINAL NUMBER OF EIGENVALUES BEING SOUGHT INTEGER ONMSGT * THE ORIGINAL SHIFT BEING SEARCHED AROUND DOUBLE PRECISION OSIGMA * THE CURRENT LEFT SHIFT DOUBLE PRECISION LFTVAL * THE CURRENT RIGHT SHIFT DOUBLE PRECISION RGTVAL * DO WE HAVE A LEFT SHIFT LOGICAL NWLEFT * DO WE HAVE A RIGHT SHIFT LOGICAL NWRGHT * THE DEBUGGING VALUE INTEGER DEBUG * THE NORMS OF K AND M DOUBLE PRECISION KNORM, MNORM * THE MACHINE EPSILON DOUBLE PRECISION EPS * THE TOLERANCE NEEDED FOR EIGENVALUE ACCEPTANCE DOUBLE PRECISION ATOL * THE ERROR CHECKING LEVEL INTEGER ERRCHK * THE CURRENT NUMBER OF EIGENVALUES BEING SOUGHT INTEGER NMSGT * DOUBLE PRECISION WORK VECTOR DOUBLE PRECISION WORKV1(*) * THE EIGENVALUES IN SORTED ORDER DOUBLE PRECISION SRTTHT(*) * THE RANK OF THE EIGENVALUE IN DISTANCE FROM OSIGMA INTEGER THTIDX(*) * A LIST OF INERTIA CALCULATIONS INTEGER INTLST(*) * A LIST OF SIGMAS DOUBLE PRECISION SIGLST(*) * THE NUMBER OF SIGMAS AND INERTIA INTEGER INTIDX * THE ACCESSING INDEX FOR FA INTEGER FALEAD * STORAGE FOR THE FACTORED MATRIX DOUBLE PRECISION FA(*) * STORAGE FOR THE FACTORED MATRIX INTEGER FARP(*) * STORAGE FOR THE FACTORED MATRIX INTEGER FARL(*) * FACTORIZATION WORK VECTORS INTEGER FIWRKV(N+1,2) DOUBLE PRECISION FRWRKV(N,3) * VIBRATION (0) OR BUCKLING (1) INTEGER PROB * THE LEVEL OF LOOP UNROLLING INTEGER NROLL * THE STIFFNESS MATRIX DOUBLE PRECISION KAUX(*) * THE STIFFNESS MATRIX INTEGER KRP(*) * THE STIFFNESS MATRIX INTEGER KCP(*) * THE MASS MATRIX DOUBLE PRECISION MAUX(*) * THE MASS MATRIX INTEGER MRP(*) * THE MASS MATRIX INTEGER MCP(*) * THE SIZE OF FA INTEGER FASIZE * THE DOUBLE PRECISION ADDRESS OF FA INTEGER FAADDR * THE TOTAL MEMORY ALLOCATED INTEGER TMEM * THE FORM THAT K IS STORED IN INTEGER KFORM * THE FORM THAT M IS STORED IN INTEGER MFORM * THE TYPE OF FACTORIZATION INTEGER FACTYP * POINTER TO A FACTORIZATION WORK AREA INTEGER WRKPTR * SINCE THIS IS A PROCEDURE NOW, THIS IS THE RETURN VALUE INTEGER RETVAL * THE STATE WE ARE IN INTEGER TSTATE * BOUNDARIES? INTEGER CBOUND * THE RIGHT AND LEFT BOUNDARIES DOUBLE PRECISION BLEFT, BRIGHT * THE B-K STORAGE FACTOR DOUBLE PRECISION BFACT * THE MAXIMUM NUMBER OF DELAYED PIVOTS THAT CAN BE STORED INTEGER MAXDLY * SIZE OF INDICES INTO FACTORED MATRIX (FOR SPARSE FACTOR) INTEGER INSIZE * OFFSET ADDRESS INTO INDICES FOR FACTORED MATRIX INTEGER INPTR * ADDRESS INTO INDICES FOR FACTORED MATRIX INTEGER INADDR * INTERNAL VARIABLES * COUNT VARIABLES INTEGER I * DO WE DO A POST CHECK? LOGICAL PSTERR * DISTANCE BETWEEN FURTHEST EIGENVALUE AND OSIGMA DOUBLE PRECISION EIGBND * HAVE ENOUGH EIGENVALUES BEEN FOUND? LOGICAL DONE * IS EVERYTHING OKAY? LOGICAL OKAY * HAS A NEW SHIFT BEEN FOUND? LOGICAL SHTFND * MIN OF OEIGNUM AND ONMSGT INTEGER OINDX * THE RETURN VALUE FROM CHKGAP INTEGER CHKSTS * THE MISSING VALUE AND ITS ERROR BOUND DOUBLE PRECISION MISVAL, MISBJ * HAVE WE FINISHED CHECKING? LOGICAL FINCHK * UNRECOVERABLE ERROR? LOGICAL BADEND * NUMBER OF MISSING VALUES INTEGER MISSED * THE MINIMUM DIST. A SHIFT CAN BE FROM AN EIGENVALUE DOUBLE PRECISION MNDIST * TWO TEMPORARY VALUES DOUBLE PRECISION T1VAL, T2VAL * PARAMETER PASSED TO ADJSHT DOUBLE PRECISION FRACT * THE OLD RIGHT SHIFT DOUBLE PRECISION ORGTVL * DID WE HAVE AN OLD RIGHT SHIFT LOGICAL ONWRGT * INTEGER FUNCTIONS INTEGER SRCHLW, CHKGAP, GBOUND * IF NO EIGENVALUES OR NOT IN STATE 0 THEN RETURN IF (((OEIGNM.EQ.0).AND.(CBOUND.NE.1)).OR. C ((TSTATE.GT.0).AND.(TSTATE.LT.5))) THEN RETVAL = 1 RETURN ENDIF * INITIALIZE SOME VALUES ONWRGT = NWRGHT OKAY = .FALSE. SHTFND = .FALSE. FINCHK = .FALSE. BADEND = .FALSE. NWRGHT = .FALSE. ORGTVL = RGTVAL IF (CBOUND.EQ.1) THEN PSTERR = .FALSE. OINDX = OEIGNM DONE = .TRUE. ELSE * DO WE DO A POST TEST USING INERTIA? IF (ERRCHK.GE.1) THEN PSTERR = .TRUE. ELSE PSTERR = .FALSE. ENDIF * COULD WE BE FINISHED? IF (OEIGNM.GE.ONMSGT) THEN DONE = .TRUE. OINDX = ONMSGT ELSE DONE = .FALSE. OINDX = OEIGNM ENDIF * SEARCH FOR UNCONVERGED NEEDED EIGENVALUES MISSED = SRCHLW(OEIGNM,OTHETA,R1LEN,R1THET,R1BJ, C R2LEN,R2THET,R2BJ,OINDX,OSIGMA,WORKV1, C MISVAL,MISBJ,ATOL) IF (MISSED.GT.0) THEN FINCHK = .TRUE. RGTVAL = MISVAL NMSGT = MAX(OEIGNM + MISSED,NMSGT) NMSGT = MIN(NMSGT,N) IF (DEBUG.GT.0) THEN PRINT *,'MISSED = ',MISSED,OEIGNM,NMSGT ENDIF * ADJUST NEW RIGHT SHIFT MNDIST = (1.0D0/(ATOL/RGTVAL)) * C ((KNORM - (MNORM*RGTVAL))*EPS) MNDIST = MAX(MNDIST,0.001D0*RGTVAL) MNDIST = MIN(MNDIST,0.01D0*RGTVAL) T1VAL = RGTVAL T2VAL = MISBJ FRACT = 0.1D0 IF (DEBUG.GT.0) THEN PRINT *,'RIGHT VAL = ',RGTVAL,MISBJ PRINT *,'MIN DIST = ',MNDIST ENDIF CALL ADJSHT(RGTVAL,T1VAL,T2VAL,MNDIST,OSIGMA, C OEIGNM,OTHETA,OLDBJ,R1LEN,R1THET,R1BJ, C R2LEN,R2THET,R2BJ,FRACT) IF (DEBUG.GT.0) THEN PRINT *,'ADJUSTED RIGHT VAL = ',RGTVAL ENDIF IF (ONWRGT.AND.(ABS(ORGTVL-MISVAL).GT. C ABS(RGTVAL-MISVAL))) THEN NWRGHT = .TRUE. ELSE IF (NWLEFT.AND.(ABS(LFTVAL-MISVAL).GT. C ABS(RGTVAL-MISVAL))) THEN NWRGHT = .TRUE. ELSE NWRGHT = .FALSE. BADEND = .TRUE. PRINT *,'ERROR: MACHINE PRECISION TOO LOW TO GET', C 'POSSIBLE EIGENVALUE AT ',MISVAL PRINT *,'REMEDY: EITHER SEARCH FOR MORE EIGENVALUES' PRINT *,'REMEDY: OR TRY A LOWER ERROR TOLERANCE' ENDIF ENDIF ENDIF IF (BADEND) THEN RETVAL = -1 RETURN ELSE IF (NWRGHT) THEN RETVAL = 2 RETURN ENDIF * SORT THE EIGENVALUES AND RELATED DATA CALL CPYVEC(OEIGNM,SRTTHT,OTHETA) CALL BSORT(OEIGNM,SRTTHT) * FRWRKV(1) CONTAINS THE SORTED DISTANCES OF EIGENVALUE FROM OSIGMA * FIWRKV(1) THE INDICES OF THE SORTED DISTANCES DO 10 I = 1, OEIGNM FRWRKV(I,1) = ABS(SRTTHT(I)-OSIGMA) FIWRKV(I,1) = I 10 CONTINUE CALL BSORT2(OEIGNM,FRWRKV(1,1),FIWRKV(1,1)) DO 20 I = 1, OEIGNM THTIDX(FIWRKV(I,1)) = I 20 CONTINUE * GIVE THE DISTANCE OF THE LAST SIGNIFICANT EIGENVALUE FROM OSIGMA EIGBND = FRWRKV(OINDX,1) * DONE WITH FRWRKV(1) AND FIWRKV(1) IF (CBOUND.EQ.1) THEN ONMSGT = GBOUND(INTIDX,SIGLST,INTLST,BLEFT,BRIGHT,0) ENDIF * THE MAIN LOOP EXECUTED UNTIL WE REACH A DECISION 100 CONTINUE IF ((OEIGNM.GE.ONMSGT).OR.(CBOUND.EQ.1)) THEN * CHECK GAPS CHKSTS = CHKGAP(INTIDX,SIGLST,INTLST, C OEIGNM,SRTTHT,FIWRKV(1,1),THTIDX,FRWRKV(1,1),WORKV1, C FIWRKV(1,2),0.001D0,ONMSGT,N,PSTERR,RGTVAL,OSIGMA,NMSGT, C CBOUND,BLEFT,BRIGHT,PROB) IF (DEBUG.GT.0) THEN PRINT *,'CHKSTS = ',CHKSTS ENDIF IF (CHKSTS.EQ.-1) THEN NWRGHT = .TRUE. FINCHK = .TRUE. ELSE IF (CHKSTS.EQ.-3) THEN CALL BOPER(N,WORKV1,WORKV1,RGTVAL,PROB,6,FACTYP, C FALEAD,FA,FARP,FARL,FIWRKV,FRWRKV,INTLST, C SIGLST,INTIDX,DEBUG,NROLL,WRKPTR, C FASIZE,FAADDR,TMEM,KFORM,KAUX,KRP,KCP, C MFORM,MAUX,MRP,MCP,BFACT,MAXDLY,INSIZE,INPTR,INADDR) ELSE IF (CHKSTS.EQ.-2) THEN BADEND = .TRUE. ELSE IF (CHKSTS.EQ.0) THEN OKAY = .TRUE. ELSE IF (CHKSTS.EQ.1) THEN FINCHK = .TRUE. ENDIF ELSE FINCHK = .TRUE. ENDIF IF (OKAY.OR.FINCHK.OR.BADEND) GOTO 110 GOTO 100 110 CONTINUE NWLEFT = .FALSE. IF (BADEND) THEN RETVAL = -1 ELSE IF (OKAY) THEN RETVAL = 0 ELSE IF (NWRGHT) THEN RETVAL = 2 ELSE IF (FINCHK) THEN RETVAL = 1 ELSE PRINT *,'ERROR: BUG IN CHKFIN' PRINT *,'REMEDY: CONTACT TESTBED ADMINISTRATOR' RETVAL = -1 ENDIF RETURN END * FINDS IF ANY EIGENVALUES WERE MISSED INTEGER FUNCTION SRCHLW(OEIGNM,OTHETA,R1LEN,R1THET,R1BJ, C R2LEN,R2THET,R2BJ,ONMSGT,OSIGMA,WORKV,MISVAL,MISBJ,ATOL) * NUMBER OF GOOD EIGS INTEGER OEIGNM * GOOD THETA'S DOUBLE PRECISION OTHETA(*) * # OF OTHER THETA'S INTEGER R1LEN * OTHER THETA'S DOUBLE PRECISION R1THET(*) * ERROR BOUNDS FOR R1'S DOUBLE PRECISION R1BJ(*) * # OF OTHER THETA'S INTEGER R2LEN * OTHER THETA'S DOUBLE PRECISION R2THET(*) * ERROR BOUNDS FOR R2'S DOUBLE PRECISION R2BJ(*) * THE ORIGINAL NUMBER OF EIGENVALUES SOUGHT INTEGER ONMSGT * THE ORIGINAL SIGMA TO SEARCH AROUND DOUBLE PRECISION OSIGMA * A WORK VECTOR DOUBLE PRECISION WORKV(*) * ONE OF THE MISSING VALUES DOUBLE PRECISION MISVAL * THE BJ OF THE MISSING VALUE DOUBLE PRECISION MISBJ * THE TOLERANCE NEEDED FOR EIGENVALUE ACCEPTANCE DOUBLE PRECISION ATOL * INTERNAL VARIABLES * COUNT VARIABLES INTEGER I, J * THE NUMBER OF EIGENVALUES GREATER THAN THE MISSED VALUE INTEGER LEVEL * THE TOTAL NUMBER OF VALUES LESS THAN SIGMA INTEGER TOTAL * HAVE WE FOUND A MISSING VAL YET? LOGICAL FIRST * THE MINIMUM DISTANCE FROM OSIGMA FOR A MISSED EIGENVALUE DOUBLE PRECISION MNDIST * THE DISTANCE FROM OSIGMA FOR EIGENVALUE DOUBLE PRECISION DIST DO 10 I = 1, OEIGNM WORKV(I) = ABS(OSIGMA-OTHETA(I)) 10 CONTINUE CALL BSORT(OEIGNM,WORKV) TOTAL = 0 LEVEL = 0 FIRST = .TRUE. DO 20 I = 1, R1LEN IF ((R1BJ(I).LT.0.0D0).AND.(R1BJ(I).GT.(-100.0D0*ATOL))) THEN DO 30 J = ONMSGT, 1, -1 DIST = ABS(OSIGMA-R1THET(I)) IF (DIST.LT.WORKV(J)) THEN LEVEL = MAX0((ONMSGT-J)+1,LEVEL) IF (J.EQ.ONMSGT) THEN TOTAL = TOTAL + 1 IF (FIRST) THEN FIRST = .FALSE. MISVAL = R1THET(I) MISBJ = R1BJ(I) MNDIST = DIST ELSE IF (DIST.LT.MNDIST) THEN MISVAL = R1THET(I) MISBJ = R1BJ(I) MNDIST = DIST ENDIF ENDIF ENDIF ELSE GOTO 40 ENDIF 30 CONTINUE 40 CONTINUE ENDIF 20 CONTINUE DO 50 I = 1, R2LEN IF ((R2BJ(I).LT.0.0D0).AND.(R2BJ(I).GT.(-100.0D0*ATOL))) THEN DO 60 J = ONMSGT, 1, -1 DIST = ABS(OSIGMA-R2THET(I)) IF (DIST.LT.WORKV(J)) THEN LEVEL = MAX0((ONMSGT-J)+1,LEVEL) IF (J.EQ.ONMSGT) THEN TOTAL = TOTAL + 1 IF (FIRST) THEN FIRST = .FALSE. MISVAL = R2THET(I) MISBJ = R2BJ(I) MNDIST = DIST ELSE IF (DIST.LT.MNDIST) THEN MISVAL = R2THET(I) MISBJ = R2BJ(I) MNDIST = DIST ENDIF ENDIF ENDIF ELSE GOTO 70 ENDIF 60 CONTINUE 70 CONTINUE ENDIF 50 CONTINUE SRCHLW = MIN(TOTAL,LEVEL) RETURN END * CHECKS GAPS BETWEEN INERTIA CHECKS INTEGER FUNCTION CHKGAP(INTIDX,SIGLST,INTLST, C NUMEIG,SRTTHT,THTACT,THTIDX,THTCLS,GAPNUM,GAPNM2,ATOL, C ONMSGT,N,PSTERR,SIGMA,OSIGMA,NMSGT,CBOUND,BLEFT,BRIGHT,PROB) * THE NUMBER OF SHIFTS WITH INERTIA CHECKS TAKEN INTEGER INTIDX * THE VALUES OF THE SHIFTS DOUBLE PRECISION SIGLST(*) * THE INERTIAS AT THE SHIFTS INTEGER INTLST(*) * THE NUMBER OF EIGENVALUES FOUND SO FAR INTEGER NUMEIG * THE SORTED EIGENVALUES DOUBLE PRECISION SRTTHT(*) * THE STATUS OF THE SORTED EIGENVALUES INTEGER THTACT(*) * THE RANK OF THE SORTED EIGENVALUES INTEGER THTIDX(*) * ARE THESE EIGENVALUES VERY CLOSE? LOGICAL THTCLS(*) * THE NUMBER OF EIGENVALUES IN EACH GAP INTEGER GAPNUM(*) * THE NUMBER OF USED EIGENVALUES IN EACH GAP INTEGER GAPNM2(*) * THE ERROR TOLERANCE ALLOWED DOUBLE PRECISION ATOL * THE NUMBER OF EIGENVALUES REQUESTED BY THE USER INTEGER ONMSGT * THE NUMBER OF EQUATIONS INTEGER N * ARE WE LOOKING FOR UNCHECKED GAPS? LOGICAL PSTERR * THE SHIFT TO WORK ON DOUBLE PRECISION SIGMA * THE ORIGINAL SIGMA DOUBLE PRECISION OSIGMA * THE NUMBER OF EIGENVALUES THAT LANZ IS SEARCHING FOR INTEGER NMSGT * BOUNDARIES? INTEGER CBOUND * THE LEFT AND RIGHT BOUNDARIES DOUBLE PRECISION BLEFT, BRIGHT * VIBRATION (0) OR BUCKLING (1) INTEGER PROB * INTERNAL EIGENVALUES * BOUNDARIES FOR MISSED VALUE DOUBLE PRECISION LEFT, RIGHT * DIRECTION TO SEARCH IN INTEGER DIR * STARTING VALUE TO SEARCH FROM DOUBLE PRECISION START * COUNT VARIABLES INTEGER I,J * ARE THERE UNACCOUNTED VALUES IN THE TOP AND BOTTOM GAPS? LOGICAL TZERO, BZERO * THE INDICES OF THE HI AND LOW USED EIGENVALUES INTEGER HINUM, LONUM * USED FOR DETERMINING BOUNDARY EIGENVALUES INTEGER THINUM, TLONUM * THE NUMBER OF GOOD EIGENVALUES ACCOUNTED FOR INTEGER NUMACT * IS THERE A GAP CONTAINING EIGENVALUES WE NEED LOGICAL SIGGAP * IS THERE A GAP CONTAINING TOO MANY EIGENVALUES LOGICAL BADGAP * TEMPORARY VERSION OF SIGGAP LOGICAL TMPGAP * THE NUMBER OF MISSING EIGENVALUES INTEGER SIGNUM * TEMPORARY VERSION OF SIGNUM INTEGER TMPNUM * THE MIDPOINT IN A GAP DOUBLE PRECISION MIDPNT * THE DISTANCE A GAP IS FROM OSIGMA DOUBLE PRECISION DIST * TEMPORARY VERSION OF DIST DOUBLE PRECISION TDIST * LOGICAL FUNCTIONS LOGICAL ECLOSE, GSHIFT * THE NUMBER OF EIGENVALUES EXPECTED BETWEEN TWO SHIFTS INTEGER NUMNGP IF ((NUMEIG.LT.ONMSGT).AND.(CBOUND.EQ.0)) THEN PRINT *,'ERROR: NOT ENOUGH EIGENVALUES TO CHECK GAPS' PRINT *,'REMEDY: CONTACT TESTBED ADMINISTRATOR' STOP ENDIF * SORT THE SHIFTS, INITIALIZE GAPNM, AND INITIALIZE THTACT CALL BSORT2(INTIDX,SIGLST,INTLST) DO 5 I = 1, INTIDX+1 GAPNUM(I) = 0 GAPNM2(I) = 0 5 CONTINUE DO 7 I = 1, NUMEIG THTACT(I) = -1 7 CONTINUE DO 9 I = 1, NUMEIG-1 THTCLS(I) = ECLOSE(ATOL,SRTTHT(I),SRTTHT(I+1)) 9 CONTINUE * CHECK OFF EIGENVALUES IN GAPS TLONUM = 0 LONUM = 0 DO 10 I = 1, NUMEIG IF ((LONUM.EQ.0).AND.(THTIDX(I).LE.ONMSGT)) THEN LONUM = I ENDIF DO 20 J = 2, INTIDX IF ((SRTTHT(I).GT.SIGLST(J-1)).AND. C (SRTTHT(I).LT.SIGLST(J))) THEN GAPNUM(J) = GAPNUM(J) + 1 IF (THTIDX(I).LE.ONMSGT) THEN GAPNM2(J) = GAPNM2(J) + 1 ENDIF THTACT(I) = 1 TLONUM = -1 GOTO 30 ENDIF 20 CONTINUE * THIS COULD BE A BOUNDARY EIGENVALUE IF (THTIDX(I).LE.ONMSGT) THEN IF (TLONUM.EQ.0) THEN THTACT(I) = 2 TLONUM = I ELSE IF ((TLONUM.GT.0).AND.THTCLS(I-1).AND. C (THTACT(I-1).EQ.2)) THEN THTACT(I) = 2 ELSE THTACT(I) = 0 ENDIF ELSE THTACT(I) = 0 ENDIF IF (INTIDX.EQ.0) THEN GAPNUM(1) = GAPNUM(1) + 1 IF (THTIDX(I).LE.ONMSGT) THEN GAPNM2(1) = GAPNM2(1) + 1 ENDIF ELSE IF (SRTTHT(I).LE.SIGLST(1)) THEN GAPNUM(1) = GAPNUM(1) + 1 THTACT(I) = THTACT(I) + 8 ELSE IF (SRTTHT(I).GT.SIGLST(INTIDX)) THEN GAPNUM(INTIDX+1) = GAPNUM(INTIDX+1) + 1 THTACT(I) = THTACT(I) + 16 ENDIF 30 CONTINUE 10 CONTINUE THINUM = 0 HINUM = 0 * TRAVERSE IN THE OTHER DIRECTION TO CLEAN UP THE DO 40 I = NUMEIG, 1, -1 IF ((HINUM.EQ.0).AND.(THTIDX(I).LE.ONMSGT)) THEN HINUM = I ENDIF IF (THTACT(I)/16.GE.1) THEN IF ((GAPNUM(INTIDX+2).EQ.(N-INTLST(INTIDX))).AND. C (PROB.EQ.0)) THEN THTACT(I) = 1 ELSE THTACT(I) = THTACT(I) - 16 ENDIF IF (THTIDX(I).LE.ONMSGT) THEN GAPNM2(INTIDX+1) = GAPNM2(INTIDX+1) + 1 ENDIF ELSE IF (THTACT(I)/8.GE.1) THEN IF ((GAPNUM(1).EQ.INTLST(1)).AND.(PROB.EQ.0)) THEN THTACT(I) = 1 ELSE THTACT(I) = THTACT(I) - 8 ENDIF IF (THTIDX(I).LE.ONMSGT) THEN GAPNM2(1) = GAPNM2(1) + 1 ENDIF ENDIF IF (THTIDX(I).LE.ONMSGT) THEN IF (THINUM.EQ.0) THEN IF (THTACT(I).EQ.1) THEN THINUM = -1 ELSE THINUM = I THTACT(I) = 3 ENDIF ELSE IF ((THINUM.GT.0).AND.THTCLS(I).AND. C (THTACT(I+1).EQ.3)) THEN THTACT(I) = 3 ENDIF ENDIF 40 CONTINUE IF (TLONUM.GT.0) THEN IF (THTACT(TLONUM).NE.2) THEN TLONUM = -1 ENDIF ENDIF TZERO = .FALSE. BZERO = .FALSE. DO 45 I = 1, NUMEIG IF ((THTACT(I).EQ.2).AND.(TLONUM.LT.1).AND.(PROB.EQ.0)) THEN THTACT(I) = 0 ENDIF IF ((THTACT(I).EQ.0).AND.(THTIDX(I).LE.ONMSGT)) THEN IF (INTIDX.EQ.0) THEN TZERO = .TRUE. BZERO = .TRUE. ELSE IF (SRTTHT(I).LT.SIGLST(1)) THEN BZERO = .TRUE. ELSE IF (SRTTHT(I).GT.SIGLST(INTIDX)) THEN TZERO = .TRUE. ENDIF ENDIF ENDIF 45 CONTINUE * MAKE SURE THAT NO GAP HAS TOO MANY EIGENVALUES * AND THAT NO GAP IS MISSING DESIRED EIGENVALUES BADGAP = .FALSE. SIGGAP = .FALSE. NUMACT = ONMSGT - GAPNM2(1) IF ((INTIDX.GT.0).AND.(GAPNUM(1).GT.INTLST(1)).AND.(PROB.EQ.0)) C THEN BADGAP = .TRUE. PRINT *,'ERROR: TOO MANY EIGENVALUES BETWEEN -INFINITY AND ', C SIGLST(1) PRINT *,'REMEDY: EXAMINE OUTPUT AND ', C 'CONTACT TESTBED ADMINISTRATOR' ENDIF IF ((INTIDX.GT.0).AND.(GAPNUM(INTIDX+1).GT.N-INTLST(INTIDX)).AND. C (PROB.EQ.0)) THEN BADGAP = .TRUE. PRINT *,'ERROR: TOO MANY EIGENVALUES BETWEEN ',SIGLST(INTIDX), C ' AND +INFINITY' PRINT *,'REMEDY: EXAMINE OUTPUT AND ', C 'CONTACT TESTBED ADMINISTRATOR' ENDIF DO 50 I = 2, INTIDX IF (PROB.EQ.0) THEN NUMNGP = INTLST(I) - INTLST(I-1) ELSE IF ((SIGLST(I).GT.0.0D0).AND.(SIGLST(I-1).GE.0.0D0)) THEN NUMNGP = INTLST(I) - INTLST(I-1) ELSE IF ((SIGLST(I).GT.0.0D0).AND.(SIGLST(I-1).LT.0.0D0)) C THEN NUMNGP = INTLST(I) + INTLST(I-1) ELSE NUMNGP = INTLST(I-1) - INTLST(I) ENDIF ENDIF IF (GAPNUM(I).GT.NUMNGP) THEN BADGAP = .TRUE. PRINT *,'ERROR: TOO MANY EIGENVALUES BETWEEN ', C SIGLST(I-1),' AND ',SIGLST(I) PRINT *,'REMEDY: EXAMINE OUTPUT AND ', C 'CONTACT TESTBED ADMINISTRATOR' ELSE IF (GAPNUM(I).LT.NUMNGP) THEN MIDPNT = (SIGLST(I)+SIGLST(I-1))/2.0D0 TMPNUM = NUMNGP - GAPNUM(I) IF ((CBOUND.EQ.1).AND.((SIGLST(I).GT.BRIGHT).OR. C (SIGLST(I-1).LT.BLEFT))) THEN TMPGAP = .FALSE. ELSE IF (CBOUND.EQ.1) THEN TMPGAP = .TRUE. ELSE IF ((NUMACT-GAPNM2(I).GT.0).AND.(NUMACT.LT.ONMSGT)) C THEN TMPGAP = .TRUE. ELSE IF ((NUMACT.LT.ONMSGT).AND.(OSIGMA.GE.MIDPNT)) THEN TMPGAP = .TRUE. ELSE IF ((NUMACT-GAPNM2(I).GT.0).AND. C (OSIGMA.LE.MIDPNT)) THEN TMPGAP = .TRUE. ELSE TMPGAP = .FALSE. ENDIF IF (TMPGAP) THEN IF ((OSIGMA.GT.SIGLST(I-1)).AND. C (OSIGMA.LT.SIGLST(I))) THEN TDIST = 0.0D0 ELSE TDIST = MIN(ABS(OSIGMA-SIGLST(I-1)), C ABS(OSIGMA-SIGLST(I))) ENDIF IF (SIGGAP) THEN IF (TDIST.LT.DIST) THEN LEFT = SIGLST(I-1) RIGHT = SIGLST(I) DIST = TDIST SIGNUM = TMPNUM ENDIF ELSE SIGGAP = .TRUE. DIST = TDIST LEFT = SIGLST(I-1) RIGHT = SIGLST(I) SIGNUM = TMPNUM ENDIF ENDIF ENDIF NUMACT = NUMACT - GAPNM2(I) 50 CONTINUE IF (BADGAP) THEN CHKGAP = -2 RETURN ELSE IF (SIGGAP) THEN CHKGAP = -1 NMSGT = MIN(NUMEIG + SIGNUM,N) START = (RIGHT+LEFT)/2.0D0 IF (ABS(OSIGMA-LEFT).LT.(ABS(OSIGMA-RIGHT))) THEN DIR = -1 ELSE DIR = 1 ENDIF IF (GSHIFT(SIGMA,LEFT,RIGHT,-1,START,ATOL,NUMEIG, C SRTTHT,THTIDX,ONMSGT,.FALSE.)) THEN CHKGAP = -1 RETURN ELSE PRINT *,'ERROR: THERE IS A MISSING EIGENVALUE, BUT THE', C ' EIGENVALUES WERE TOO CLUSTERED TO CONTINUE' PRINT *,'REMEDY: TRY A PUTTING A SHIFT IN THE BAD REGION', C ' YOURSELF' CHKGAP = -2 RETURN ENDIF ENDIF * CHECK FOR VALUES THAT MAY NOT BE INERTIA CHECKED IF (PSTERR) THEN IF (INTIDX.EQ.0) THEN IF (HINUM.LE.2) THEN PRINT *,'ERROR: NO PLACE FOR INERTIA CHECK SHIFT' PRINT *,'REMEDY: NEED TO FIND MORE EIGENVALUES' CHKGAP = -2 RETURN ELSE START = (SRTTHT(HINUM)+SRTTHT(HINUM-1))/2 LEFT = SRTTHT(LONUM) RIGHT = SRTTHT(HINUM) IF (GSHIFT(SIGMA,LEFT,RIGHT,-1,START,ATOL,NUMEIG, C SRTTHT,THTIDX,ONMSGT,.TRUE.)) THEN CHKGAP = -3 RETURN ELSE PRINT *,'ERROR: THE EIGENVALUES WERE TOO', C ' CLUSTERED TO DO AN INERTIA CHECK' PRINT *,'REMEDY: EITHER FIND MORE EIGENVALUES OR', C ' RESTART AND CHOOSE A NEW SHIFT YOURSELF' CHKGAP = -2 RETURN ENDIF ENDIF ELSE IF (TZERO) THEN LEFT = SIGLST(INTIDX) RIGHT = SRTTHT(HINUM) START = (SRTTHT(HINUM)+SRTTHT(HINUM-1))/2 IF (GSHIFT(SIGMA,LEFT,RIGHT,-1,START,ATOL,NUMEIG, C SRTTHT,THTIDX,ONMSGT,.TRUE.)) THEN CHKGAP = -3 RETURN ELSE PRINT *,'ERROR: THE EIGENVALUES WERE TOO', C ' CLUSTERED TO DO AN INERTIA CHECK' PRINT *,'REMEDY: EITHER FIND MORE EIGENVALUES OR', C ' RESTART AND CHOOSE A NEW SHIFT YOURSELF' CHKGAP = -2 RETURN ENDIF ELSE IF (BZERO) THEN LEFT = SRTTHT(LONUM) RIGHT = SIGLST(1) START = (SRTTHT(LONUM)+SRTTHT(LONUM+1))/2 IF (GSHIFT(SIGMA,LEFT,RIGHT,1,START,ATOL,NUMEIG, C SRTTHT,THTIDX,ONMSGT,.TRUE.)) THEN CHKGAP = -3 RETURN ELSE PRINT *,'ERROR: THE EIGENVALUES WERE TOO', C ' CLUSTERED TO DO AN INERTIA CHECK' PRINT *,'REMEDY: EITHER FIND MORE EIGENVALUES OR', C ' RESTART AND CHOOSE A NEW SHIFT YOURSELF' CHKGAP = -2 RETURN ENDIF ELSE DO 60 I = 2, INTIDX IF (PROB .EQ. 0) THEN NUMNGP = INTLST(I) - INTLST(I-1) ELSE IF ((SIGLST(I).GT.0.0D0).AND.(SIGLST(I-1).GE.0.0D0)) C THEN NUMNGP = INTLST(I) - INTLST(I-1) ELSE IF ((SIGLST(I).GT.0.0D0).AND. C (SIGLST(I-1).LT.0.0D0)) THEN NUMNGP = INTLST(I) + INTLST(I-1) ELSE NUMNGP = INTLST(I-1) - INTLST(I) ENDIF ENDIF IF ((GAPNUM(I).LT.NUMNGP).AND. C (GAPNM2(I).GT.0)) THEN IF (GAPNM2(I).EQ.1) THEN IF ((SRTTHT(HINUM).GT.SIGLST(I-1)).AND. C(SRTTHT(HINUM).LT.SIGLST(I))) THEN GOTO 60 ENDIF IF ((SRTTHT(LONUM).GT.SIGLST(I-1)).AND. C(SRTTHT(LONUM).LT.SIGLST(I))) THEN GOTO 60 ENDIF ENDIF IF ((SRTTHT(HINUM).GT.SIGLST(I-1)).AND. C(SRTTHT(HINUM).LE.SIGLST(I))) THEN IF (HINUM.EQ.1) THEN START = (SRTTHT(HINUM)+SIGLST(I-1))/2.0D0 ELSE START = MAX(SIGLST(I-1),SRTTHT(HINUM-1)) START = (SRTTHT(HINUM)+START)/2.0D0 ENDIF RIGHT = SRTTHT(HINUM) LEFT = SIGLST(I-1) DIR = -1 ELSE IF ((SRTTHT(LONUM).GE.SIGLST(I-1)).AND. C (SRTTHT(LONUM).LT.SIGLST(I))) THEN IF (LONUM.EQ.NUMEIG) THEN START = (SRTTHT(LONUM)+SIGLST(I))/2.0D0 ELSE START = MIN(SIGLST(I),SRTTHT(LONUM+1)) START = (SRTTHT(LONUM)+START)/2.0D0 ENDIF RIGHT = SIGLST(I) LEFT = SRTTHT(LONUM) DIR = 1 ELSE RIGHT = SIGLST(I) LEFT = SIGLST(I-1) START = (RIGHT+LEFT)/2.0D0 IF (ABS(OSIGMA-LEFT).LT.(ABS(OSIGMA-RIGHT))) THEN DIR = -1 ELSE DIR = 1 ENDIF ENDIF IF (GSHIFT(SIGMA,LEFT,RIGHT,DIR,START,ATOL,NUMEIG, C SRTTHT,THTIDX,ONMSGT,.TRUE.)) THEN CHKGAP = -3 RETURN ELSE PRINT *,'ERROR: THE EIGENVALUES WERE TOO', C ' CLUSTERED TO DO AN INERTIA CHECK' PRINT *,'REMEDY: EITHER FIND MORE EIGENVALUES OR', C ' RESTART AND CHOOSE A NEW SHIFT YOURSELF' CHKGAP = -2 RETURN ENDIF ENDIF 60 CONTINUE ENDIF ELSE CHKGAP = 1 ENDIF CHKGAP = 0 RETURN END LOGICAL FUNCTION ECLOSE(ATOL,EIG1,EIG2) * THE ACCURACY OF THE EIGENVALUES DOUBLE PRECISION ATOL * THE TWO EIGENVALUES DOUBLE PRECISION EIG1, EIG2 * INTERNAL VARIABLES * THE ERROR IN EACH EIGENVALUE DOUBLE PRECISION ERR1, ERR2 ERR1 = ATOL*EIG1 ERR2 = ATOL*EIG2 IF (ERR1+ERR2.GT.ABS(EIG2-EIG1)) THEN ECLOSE = .TRUE. ELSE ECLOSE = .FALSE. ENDIF RETURN END LOGICAL FUNCTION GSHIFT(SIGMA,LEFT,RIGHT,DIR,START,ATOL, C NUMEIG,SRTTHT,THTIDX,ONMSGT,GUESS) * THE NEW SHIFT DOUBLE PRECISION SIGMA * THE LEFT AND RIGHT BOUNDARIES DOUBLE PRECISION LEFT, RIGHT * THE DIRECTIONS PREFERRED INTEGER DIR * STARTING PLACE TO LOOK FOR THE SHIFT DOUBLE PRECISION START * THE ERROR TOLERANCE DOUBLE PRECISION ATOL * THE NUMBER OF EIGENVALUES FOUND INTEGER NUMEIG * THE EIGENVALUES IN SORTED ORDER DOUBLE PRECISION SRTTHT(*) * THE INDICES OF THE EIGENVALUES INTEGER THTIDX(*) * THE NUMBER OF EIGENVALUES REQUESTED BY THE USER INTEGER ONMSGT * DO WE WANT TO PUT SIGMA NEAR A BOUNDARY LOGICAL GUESS * INTERNAL VARIABLES * COUNT VARIABLES INTEGER I * THE INDEX OF THE EIGENVALUES JUST ABOVE START INTEGER MINDEX * THE VALUE OF THE MINDEX'TH EIGENVALUES DOUBLE PRECISION MVAL * THE CURRENT EIGENVALUE BEING LOOKED AT INTEGER CINDEX * LOGICAL FUNCTIONS LOGICAL SOKAY IF (NUMEIG.EQ.0) THEN IF (SOKAY(SIGMA,LEFT,SRTTHT(CINDEX), C DIR,LEFT,RIGHT,ATOL,GUESS)) THEN GSHIFT = .TRUE. ELSE GSHIFT = .FALSE. ENDIF RETURN ENDIF GSHIFT = .FALSE. * GET BOUND MINDEX = 0 DO 10 I = 1, NUMEIG IF (SRTTHT(I).GT.START) THEN IF (MINDEX.EQ.0) THEN MINDEX = I MVAL = SRTTHT(I) ELSE IF (SRTTHT(I).LT.MVAL) THEN MINDEX = I MVAL = SRTTHT(I) ENDIF ENDIF 10 CONTINUE IF (MINDEX.EQ.0) THEN MINDEX = NUMEIG + 1 ENDIF * DECIDE WHICH TO DO FIRST IF (DIR.EQ.-1) GOTO 200 * DO RIGHT 100 CONTINUE CINDEX = MINDEX 110 CONTINUE IF (CINDEX.LT.1) GOTO 190 IF (CINDEX.GT.NUMEIG+1) GOTO 190 IF (THTIDX(CINDEX).GT.ONMSGT) GOTO 190 IF (CINDEX.EQ.1) THEN IF (SOKAY(SIGMA,LEFT,SRTTHT(CINDEX), C DIR,LEFT,RIGHT,ATOL,GUESS)) THEN GSHIFT = .TRUE. RETURN ENDIF ELSE IF (CINDEX.EQ.NUMEIG+1) THEN IF (SOKAY(SIGMA,SRTTHT(CINDEX-1),RIGHT, C DIR,LEFT,RIGHT,ATOL,GUESS)) THEN GSHIFT = .TRUE. RETURN ENDIF ELSE IF (SOKAY(SIGMA,SRTTHT(CINDEX-1),SRTTHT(CINDEX), C DIR,LEFT,RIGHT,ATOL,GUESS)) THEN GSHIFT = .TRUE. RETURN ENDIF ENDIF CINDEX = CINDEX - 1 GOTO 110 * QUIT IF LEFT FIRST ELSE DO LEFT 190 CONTINUE IF (DIR.EQ.-1) RETURN * DO LEFT 200 CONTINUE CINDEX = MINDEX 210 CONTINUE IF (CINDEX.LT.1) GOTO 290 IF (CINDEX.GT.NUMEIG+1) GOTO 290 IF (THTIDX(CINDEX).GT.ONMSGT) GOTO 290 IF (CINDEX.EQ.1) THEN IF (SOKAY(SIGMA,LEFT,SRTTHT(CINDEX), C DIR,LEFT,RIGHT,ATOL,GUESS)) THEN GSHIFT = .TRUE. RETURN ENDIF ELSE IF (CINDEX.EQ.NUMEIG+1) THEN IF (SOKAY(SIGMA,SRTTHT(CINDEX-1),RIGHT, C DIR,LEFT,RIGHT,ATOL,GUESS)) THEN GSHIFT = .TRUE. RETURN ENDIF ELSE IF (SOKAY(SIGMA,SRTTHT(CINDEX-1),SRTTHT(CINDEX), C DIR,LEFT,RIGHT,ATOL,GUESS)) THEN GSHIFT = .TRUE. RETURN ENDIF ENDIF CINDEX = CINDEX + 1 GOTO 210 * QUIT IF RIGHT FIRST ELSE DO RIGHT 290 CONTINUE IF (DIR.EQ.-1) GOTO 100 RETURN END LOGICAL FUNCTION SOKAY(SIGMA,LEFT,RIGHT,DIR,BLEFT,BRIGHT,ATOL, C GUESS) * THE NEW SHIFT DOUBLE PRECISION SIGMA * THE LEFT AND RIGHT EIGENVALUES DOUBLE PRECISION LEFT, RIGHT * THE PREFERRED DIRECTION INTEGER DIR * THE LEFT AND RIGHT BOUNDARIES DOUBLE PRECISION BLEFT, BRIGHT * THE ERROR TOLERANCE DOUBLE PRECISION ATOL * DO WE WANT TO PUT SIGMA NEAR A BOUNDARY LOGICAL GUESS * INTERNAL VARIABLES * THE BOUNDARIES ON SIGMA DOUBLE PRECISION TLEFT, TRIGHT * TOLERANCE ON THE BOUNDS DOUBLE PRECISION BTOL BTOL = 0.0001D0 SOKAY = .FALSE. TLEFT = MAX((LEFT + BTOL*ABS(LEFT)),(BLEFT + BTOL*ABS(BLEFT))) TRIGHT = MIN((RIGHT - BTOL*ABS(RIGHT)), C (BRIGHT - BTOL*ABS(BRIGHT))) IF (TLEFT.LT.TRIGHT) THEN IF (GUESS) THEN IF (DIR.EQ.-1) THEN SIGMA = TRIGHT SOKAY = .TRUE. ELSE SIGMA = TLEFT SOKAY = .TRUE. ENDIF ELSE SIGMA = (TLEFT+TRIGHT)/2.0D0 SOKAY = .TRUE. ENDIF ENDIF RETURN END C$FORTRAN DBLAS *********************************************************************** * LANZ SOFTWARE PACKAGE * * FILENAME: DBLAS.MSC * * AUTHOR: RECEIVED FROM NETLIB * * PURPOSE: PERFORM VECTOR-VECTOR OPERATIONS * *********************************************************************** DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) C C TAKES THE SUM OF THE ABSOLUTE VALUES. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DTEMP INTEGER I,INCX,M,MP1,N,NINCX C DASUM = 0.0D0 DTEMP = 0.0D0 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 DTEMP = DTEMP + DABS(DX(I)) 10 CONTINUE DASUM = DTEMP RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,6) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DABS(DX(I)) 30 CONTINUE IF( N .LT. 6 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,6 DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I + 1)) + DABS(DX(I + 2)) * + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5)) 50 CONTINUE 60 DASUM = DTEMP RETURN END SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DA INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF (DA .EQ. 0.0D0) 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 DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1) DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 50 CONTINUE RETURN END SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1) INTEGER I,INCX,INCY,IX,IY,M,MP1,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 DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOT = 0.0D0 DTEMP = 0.0D0 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 DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT = DTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE 60 DDOT = DTEMP RETURN END DOUBLE PRECISION FUNCTION DMACH(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 DOUBLE PRECISION EPS,TINY,HUGE,S C EPS = 1.0D0 10 EPS = EPS/2.0D0 S = 1.0D0 + EPS IF (S .GT. 1.0D0) GO TO 10 EPS = 2.0D0*EPS C S = 1.0D0 20 TINY = S S = S/16.0D0 IF (S*1.0 .NE. 0.0D0) GO TO 20 TINY = (TINY/EPS)*100.0 HUGE = 1.0D0/TINY C IF (JOB .EQ. 1) DMACH = EPS IF (JOB .EQ. 2) DMACH = TINY IF (JOB .EQ. 3) DMACH = HUGE RETURN END DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) INTEGER NEXT DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE DATA ZERO, ONE /0.0D0, 1.0D0/ C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() 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 DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(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() DOUBLE PRECISION 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 / 8.232D-11, 1.304D19 / C IF(N .GT. 0) GO TO 10 DNRM2 = 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( DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF( DX(I) .EQ. ZERO) GO TO 200 IF( DABS(DX(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 / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF( DABS(DX(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( DABS(DX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 C 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 C C C PREPARE FOR PHASE 3. C 75 SUM = (SUM * XMAX) * XMAX C C C FOR DOUBLE PRECISION 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(DABS(DX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 DNRM2 = DSQRT( 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 DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END SUBROUTINE DROT (N,DX,INCX,DY,INCY,C,S) C C APPLIES A PLANE ROTATION. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP,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 DTEMP = C*DX(IX) + S*DY(IY) DY(IY) = C*DY(IY) - S*DX(IX) DX(IX) = DTEMP 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 DTEMP = C*DX(I) + S*DY(I) DY(I) = C*DY(I) - S*DX(I) DX(I) = DTEMP 30 CONTINUE RETURN END SUBROUTINE DROTG(DA,DB,C,S) C C CONSTRUCT GIVENS PLANE ROTATION. C JACK DONGARRA, LINPACK, 3/11/78. C MODIFIED 9/27/86. C DOUBLE PRECISION DA,DB,C,S,ROE,SCALE,R,Z C ROE = DB IF( DABS(DA) .GT. DABS(DB) ) ROE = DA SCALE = DABS(DA) + DABS(DB) IF( SCALE .NE. 0.0D0 ) GO TO 10 C = 1.0D0 S = 0.0D0 R = 0.0D0 GO TO 20 10 R = SCALE*DSQRT((DA/SCALE)**2 + (DB/SCALE)**2) R = DSIGN(1.0D0,ROE)*R C = DA/R S = DB/R 20 Z = S IF( DABS(C) .GT. 0.0D0 .AND. DABS(C) .LE. S ) Z = 1.0D0/C DA = R DB = Z RETURN END SUBROUTINE DSCAL(N,DA,DX,INCX) C C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DA,DX(1) INTEGER I,INCX,M,MP1,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 DX(I) = DA*DX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END SUBROUTINE DSWAP (N,DX,INCX,DY,INCY) C C INTERCHANGES TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,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 DTEMP = DX(IX) DX(IX) = DY(IY) DY(IY) = DTEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,3) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DX(I) DX(I) = DY(I) DY(I) = DTEMP 30 CONTINUE IF( N .LT. 3 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,3 DTEMP = DX(I) DX(I) = DY(I) DY(I) = DTEMP DTEMP = DX(I + 1) DX(I + 1) = DY(I + 1) DY(I + 1) = DTEMP DTEMP = DX(I + 2) DX(I + 2) = DY(I + 2) DY(I + 2) = DTEMP 50 CONTINUE RETURN END INTEGER FUNCTION IDAMAX(N,DX,INCX) C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DMAX INTEGER I,INCX,IX,N C IDAMAX = 0 IF( N .LT. 1 ) RETURN IDAMAX = 1 IF(N.EQ.1)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 DMAX = DABS(DX(1)) IX = IX + INCX DO 10 I = 2,N IF(DABS(DX(IX)).LE.DMAX) GO TO 5 IDAMAX = I DMAX = DABS(DX(IX)) 5 IX = IX + INCX 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DMAX = DABS(DX(1)) DO 30 I = 2,N IF(DABS(DX(I)).LE.DMAX) GO TO 30 IDAMAX = I DMAX = DABS(DX(I)) 30 CONTINUE RETURN END C$FORTRAN DPREP *********************************************************************** * LANZ SOFTWARE PACKAGE (TEST PROGRAM) * * FILENAME: DPREP.MSC * * AUTHOR: MARK JONES * * PURPOSE: READS IN DATA FOR THE TEST PROGRAM * *********************************************************************** SUBROUTINE DPREP(IPARAM,RPARAM,IOPAR,LIST,SLIST,ALIST,TMEM, C KAUX,KRP,KCP,MAUX,MRP,MCP,RTOJ) * THESE PARAMETERS ARE BETTER DETAILED IN MAIN.MSC * IPARAM IS A SET OF INTEGER PARAMETERS FOR LANZ INTEGER IPARAM(*) * RPARAM IS A SET OF DOUBLE PRECISION PARAMETERS FOR LANZ DOUBLE PRECISION RPARAM(*) * IOPAR IS THE SET OF I/O RELATED PARAMETERS INTEGER IOPAR(*) * THE LIST OF MEMORY LOCATIONS INTEGER LIST(*) * THE LIST OF MEMORY SIZES INTEGER SLIST(*) * THE LIST OF MEMORY ADDRESSES INTEGER ALIST(*) * THE TOTAL MEMORY ALLOCATED SO FAR INTEGER TMEM * THE K MATRIX DOUBLE PRECISION KAUX(*) * THE ROW POINTERS INTEGER KRP(*) * THE COLUMN INDICES INTEGER KCP(*) * THE M MATRIX DOUBLE PRECISION MAUX(*) * THE ROW POINTERS INTEGER MRP(*) * THE COLUMN INDICES INTEGER MCP(*) * THE ROW TO JOINT VECTOR INTEGER RTOJ(*) * INTERNAL VARIABLES * A LOOP COUNT VARIABLE INTEGER I * LIBRARY AND DATASET SEQUENCE NUMBER (IGNORE IF NOT IN TESTBED) INTEGER LDM, MSDN * ARRAYS OF INPUT PARAMETERS TO BE READ IN INTEGER IV(3,30) DOUBLE PRECISION RV(3,30) * READ IN INPUT PARAMETERS FROM FORT.9 * THE PARAMETERS ARE EXPLAINED IN DRIVER.DOC * PROB READ(9,*) IV(3,7) * MXIT READ(9,*) IV(3,8) * MXLI READ(9,*) IV(3,9) * YSTO READ(9,*) IV(3,10) * CONV READ(9,*) RV(3,12) * DBUG READ(9,*) IV(3,13) * SHIF READ(9,*) RV(3,14) * PLVL READ(9,*) IV(3,15) * CHEC READ(9,*) IV(3,16) * NREQ READ(9,*) IV(3,18) * MCAS READ(9,*) IV(3,19) * FACT READ(9,*) IV(3,22) * NOSH READ(9,*) IV(3,23) * WRIT READ(9,*) IV(3,24) * READ READ(9,*) IV(3,25) * SEAR READ(9,*) IV(3,26) * LEFT READ(9,*) RV(3,27) * RIGH READ(9,*) RV(3,28) * BSTO READ(9,*) RV(3,29) * NUMD READ(9,*) IV(3,30) * SET MAXIMUM NUMBER OF EIGENVALUES THAT CAN BE STORED (YSTO) IPARAM(2) = IV(3,10) * SET THE NUMBER OF EIGENVALUES BEGIN SOUGHT (NREQ) IPARAM(3) = IV(3,18) * SET THE MAXIMUM TOTAL NUMBER OF STEPS TO TAKE (MXIT) IPARAM(4) = IV(3,8) * SET THE LEVEL OF DEBUGGING (DBUG) IPARAM(7) = IV(3,13) * SET THE TYPE OF PROBLEM (PROB) IPARAM(8) = IV(3,7) * SET THE LEVEL OF ERROR CHECKING (PLVL AND CHEC) IF (IV(3,15).EQ.0) THEN * LEVEL 0 = PRINT NOTHING IPARAM(10) = 0 ELSE IF (IV(3,15).EQ.1) THEN * LEVEL 1 = PRINT EIGENVALUES AND INERTIAS IPARAM(10) = 64 + 4 ELSE IF (IV(3,15).EQ.2) THEN * LEVEL 2 = PRINT EIGENVALUES, ESTIMATED ERRORS AND INERTIAS IPARAM(10) = 32 + 4 ELSE IF (IV(3,15).EQ.3) THEN * LEVEL 3 = PRINT EIGENVALUES, ESTIMATED ERRORS, * CALCULATED ERROR AND INERTIAS IPARAM(10) = 16 + 4 ELSE IF (IV(3,15).EQ.4) THEN * LEVEL 4 = PRINT EIGENVALUES, ESTIMATED ERRORS, * CALCULATED ERROR, Y-ORTHOGONALITY AND INERTIAS IPARAM(10) = 16 + 4 + 2 ELSE IF (IV(3,15).EQ.5) THEN * LEVEL 5 = PRINT EIGENVALUES, ESTIMATED ERRORS, * CALCULATED ERROR, Y-ORTHOGONALITY AND INERTIAS * AND THROW IN SOME FREQUENCIES IPARAM(10) = 128 + 16 + 4 + 2 ENDIF * DO AN INERTIA CHECK? IF (IV(3,16).GT.0) THEN IPARAM(9) = 1 ELSE IPARAM(9) = 0 ENDIF * SET THE MAXIMUM TOTAL NUMBER OF STEPS TO TAKE ON ONE LANCZOS * RUN (MXLI) IF (IV(3,9).EQ.0) THEN IF (IPARAM(8).EQ.0) THEN IPARAM(11) = 35 ELSE IPARAM(11) = 45 ENDIF ELSE IPARAM(11) = IV(3,9) ENDIF * SET THE WAY TO READ AND STORE THE MATRICES IPARAM(13) = 0 IPARAM(14) = 0 * SET THE LEVEL OF LOOP UNROLLING (LOOP)(ALWAYS 6 NOW) IPARAM(15) = 6 * SET THE CONSTRAINT CASE NUMBER (IGNORE IF NOT IN TESTBED) IOPAR(2) = IV(3,20) * SET THE SET NUMBER (IGNORE IF NOT IN TESTBED) IOPAR(3) = IV(3,21) * SET THE OUTPUT LIBRARY NUMBER (IGNORE IF NOT IN TESTBED) IOPAR(4) = IV(3,17) * SET THE TYPE OF FACTORIZATION IPARAM(16) = IV(3,22) * SET WHETHER DYNAMIC SHIFTING IS TURNED OFF IPARAM(17) = IV(3,23) * SET TO NO INITIAL GUESS IPARAM(18) = 0 * SET WHETHER WE ARE SEARCHING IN BOUNDARIES OR NOT IPARAM(12) = IV(3,26) * SET THE SHIFT (SHIF) RPARAM(1) = RV(3,14) * SET THE RELATIVE ERROR VALUE (CONV) RPARAM(2) = RV(3,12) * SET THE LEFT BOUNDARY RPARAM(3) = RV(3,27) * SET THE RIGHT BOUNDARY RPARAM(4) = RV(3,28) * SET THE B-K STORAGE FACTOR RPARAM(5) = RV(3,29) * SET THE NUMBER OF DELAYED PIVOTS THAT CAN BE STORED IPARAM(20) = IV(3,30) * SET THE NUMBER OF Y VECTORS TO STORE IF NONE HAS BEEN SET IF (IPARAM(2).LE.0) THEN IPARAM(2) = 1.5 * IPARAM(3) IPARAM(2) = MAX0(IPARAM(2),10) ENDIF IF (IPARAM(3).GT.IPARAM(2)) THEN PRINT *,'WARNING: MAX EIGS SOUGHT TOO LARGE' IPARAM(3) = IPARAM(2) PRINT *,'FIXED: NUMBER RESET TO ',IPARAM(2) ENDIF * READ IN SPARSE MATRICES * NOW LET'S TRY TO READ IN THE DATASETS THAT WE NEED * FIRST GET THE 'K' MATRIX * CALL RDGSP (FROM READIT.MSC) TO READ IN THE K MATRIX CALL RDGSP(LDM,MSDN,KAUX,KRP,KCP,TMEM,LIST(6),LIST(7), C LIST(8),RTOJ,LIST(12),.FALSE.,IPARAM(1),IOPAR(1),IOPAR(5), C IV(3,24),2,IV(3,25),7,SLIST(1),SLIST(2),SLIST(3), C ALIST(6),ALIST(7),ALIST(8),ALIST(12)) * IF WE ARE WORKING WITH A DIAGONAL M MATRIX THEN READ IT IN IF ((IPARAM(8).EQ.0).AND.(IV(3,19).EQ.0)) THEN * GET THE 'DEM' MATRIX * ALLOCATE SPACE FOR THE NONZEROES OF M SLIST(4) = 2*(IPARAM(1)+1) CALL FALLOC(SLIST(4),0,MAUX,TMEM,LIST(9),ALIST(9)) * IF WE ARE READING FROM THE TESTBED IF (IV(3,25).EQ.0) THEN * IF WE ARE READING IN ASCII FORMAT FROM FORT.8 ELSE IF (IV(3,25).EQ.1) THEN * READ IN THE DIAGONAL DO 67 I = 1, IPARAM(1) READ(8,*) MAUX(LIST(9)+I-1) 67 CONTINUE * IF WE ARE DOING AN UNFORMATTED READ FROM FORT.8 ELSE IF (IV(3,25).EQ.2) THEN READ(8) (MAUX(LIST(9)+I-1),I=1,IPARAM(1)) ENDIF * ALLOCATE SPACE FOR THE ROW POINTERS SLIST(5) = 2*(IPARAM(1)+1) CALL FALLOC(SLIST(5),1,MRP,TMEM,LIST(10),ALIST(10)) * ALLOCATE SPACE FOR THE COLUMN NUMBERS SLIST(6) = 2*IPARAM(1) CALL FALLOC(SLIST(6),1,MCP,TMEM,LIST(11),ALIST(11)) * SET ALL THE ROW POINTERS TO ONE DO 65 I = 1, IPARAM(1)+1 MRP(I+LIST(10)-1) = 1 65 CONTINUE * INITIALIZE THE ONLY COLUMN NUMBER TO 1 MCP(LIST(11)) = 1 IF (IV(3,24).EQ.1) THEN DO 66 I = 1, IPARAM(1) WRITE(3,*) MAUX(LIST(9)+I-1) 66 CONTINUE ELSE IF (IV(3,24).EQ.2) THEN WRITE(3) (MAUX(LIST(9)+I-1),I=1,IPARAM(1)) ENDIF * IF WE ARE USING A CONSISTENT MASS MATRIX, A SPARSE MATRIX * THEN READ IT IN ELSE IF (IPARAM(8).EQ.0) THEN * GET THE 'CEM' MATRIX * CALL RDGSP (FROM READIT.MSC) TO READ IN THE M MATRIX CALL RDGSP(LDM,MSDN,MAUX,MRP,MCP,TMEM,LIST(9),LIST(10), C LIST(11),RTOJ,LIST(12),.TRUE.,IPARAM(1),IOPAR(1),IOPAR(5), C IV(3,24),3,IV(3,25),8,SLIST(4),SLIST(5),SLIST(6), C ALIST(9),ALIST(10),ALIST(11),ALIST(12)) * IF WE ARE USING A POSSIBLY INDEFINITE MATRIX IN THE RHS OF THE * EIGENVALUE PROBLEM (SUCH AS IN THE BUCKLING PROBLEM) * THEN READ IT IN (OFTEN THIS MATRIX IS CALLED THE KG MATRIX) ELSE * GET THE 'KG' MATRIX * CALL RDGSP (FROM READIT.MSC) TO READ IN THE MATRIX CALL RDGSP(LDM,MSDN,MAUX,MRP,MCP,TMEM,LIST(9),LIST(10), C LIST(11),RTOJ,LIST(12),.TRUE.,IPARAM(1),IOPAR(1),IOPAR(5), C IV(3,24),3,IV(3,25),8,SLIST(4),SLIST(5),SLIST(6), C ALIST(9),ALIST(10),ALIST(11),ALIST(12)) ENDIF RETURN END C$FORTRAN MTJEIG *********************************************************************** * LANZ SOFTWARE PACKAGE * * FILENAME: EIGCLC.MSC * * AUTHOR: MARK JONES * * PURPOSE: A TRIDIAGONAL EIGENSOLVER CALLED BY LANCZS THAT IS * * EXPLAINED IN JONES-PATRICK * *********************************************************************** SUBROUTINE EIGCLC(J,ALPHA,BET2,THETA,BJ,UNDONE, C ANORM,EPS,DEBUG,HUGE,NTHETA,NBJ) * THE SIZE OF T(J) INTEGER J * ALPHA IS THE DIAGONAL OF T(J) DOUBLE PRECISION ALPHA(*) * BET2 IS THE SQUARE OF THE OFF-DIAGONAL OF T(J) DOUBLE PRECISION BET2(*) * THETA ARE THE RITZ VALUES DOUBLE PRECISION THETA(*) * THE ERROR BOUNDS ON THETA DOUBLE PRECISION BJ(*) * AN ARRAY FOR USE IN THIS ROUTINE ONLY INTEGER UNDONE(*) * THE NORM OF THE MATRIX PENCIL DOUBLE PRECISION ANORM * EPS IS THE MACHINE PRECISION DOUBLE PRECISION EPS * DEBUG LEVEL INTEGER DEBUG * THE LARGEST DOUBLE PRECISION NUMBER DOUBLE PRECISION HUGE * THE NEW RITZ VALUES DOUBLE PRECISION NTHETA(*) * THE NEW ERROR BOUNDS DOUBLE PRECISION NBJ(*) * THE FOLLOWING ARE INTERNAL VARIABLES * THE PROBE FOR NUMLES DOUBLE PRECISION PROBE * HOLDS VALUE FROM NUMHERE INTEGER NMHERE * IS THE EIGENVALUE ISOLATED? LOGICAL OKAY * MINIMUM ERROR TOLERANCE DOUBLE PRECISION ERRTOL * TEMPORARILY HOLDS BJ VALUE DOUBLE PRECISION TEMPBJ * COUNT VARIABLE INTEGER I * ENDS OF INTERVAL DOUBLE PRECISION LEFT, RIGHT * SQRT OF BET2(J+1) DOUBLE PRECISION BETP1 * INTEGER FUNCTIONS INTEGER NUMLES EXTERNAL NUMLES, NEWTON BETP1 = SQRT(BET2(J+1)) * RETURNS EIGENVALUES AND BOUND FROM J=1 OR 2 IF (J.EQ.1) THEN THETA(1) = ALPHA(1) BJ(1) = 1.0D0*BETP1 ELSE IF (J.EQ.2) THEN THETA(1) = (ALPHA(1) + ALPHA(2) - SQRT(4.0D0*BET2(2)+ C (ALPHA(1)-ALPHA(2))*(ALPHA(1)-ALPHA(2))))/2.0D0 THETA(2) = ALPHA(1) + ALPHA(2) - THETA(1) BJ(1) = 1.0D0/(1.0D0+BET2(2)/ C (THETA(1)-ALPHA(1))*(THETA(1)-ALPHA(1))) BJ(2) = 1.0D0/(1.0D0+BET2(2)/ C (THETA(2)-ALPHA(1))*(THETA(2)-ALPHA(1))) BJ(1) = BJ(1)*BETP1 BJ(2) = BJ(2)*BETP1 ENDIF IF (J.GT.2) THEN DO 10 I = 1, J UNDONE(I) = 0 10 CONTINUE * FIND ISOLATED EIGENVALUES AND THE SPOT IN WHICH THEY BELONG DO 20 I = 1, J-1 ERRTOL = ABS(ANORM*SQRT(EPS)*THETA(I)) IF (DEBUG.GT.1) THEN PRINT *,'PRETHET = ',THETA(I),BJ(I),ERRTOL ENDIF IF (BJ(I).LT.ERRTOL) THEN TEMPBJ = ERRTOL ELSE TEMPBJ = BJ(I) ENDIF TEMPBJ = MAX(TEMPBJ + 128.0D0*TEMPBJ*EPS, C ABS(THETA(I)*(16.0D0*EPS))) IF (DEBUG.GT.1) THEN PRINT *,TEMPBJ ENDIF OKAY = .TRUE. IF (I.LT.J-1) THEN IF (2.0D0*TEMPBJ.GE.ABS(THETA(I)-THETA(I+1))) THEN OKAY = .FALSE. ENDIF ENDIF IF (I.GT.1) THEN IF (2.0D0*TEMPBJ.GT.ABS(THETA(I)-THETA(I-1))) THEN OKAY = .FALSE. ENDIF ENDIF IF (OKAY) THEN PROBE = THETA(I) + TEMPBJ NMHERE = NUMLES(ALPHA,BET2,PROBE,J,1,EPS) IF (DEBUG.GT.1) THEN PRINT *,'NUMLES = ',NMHERE,' AT ',I,PROBE ENDIF IF (NMHERE.EQ.I) THEN IF (UNDONE(I).GT.0) THEN IF (DEBUG.GT.1) THEN PRINT *,'OVERWRITING IN T(J) FINDER ',I,PROBE ENDIF ELSE UNDONE(I) = I ENDIF ELSE IF (NMHERE.EQ.I+1) THEN IF (UNDONE(I+1).GT.0) THEN IF (DEBUG.GT.1) THEN PRINT *,'OVERWRITING IN T(J) FINDER ',I,PROBE ENDIF ELSE UNDONE(I+1) = I ENDIF ELSE PRINT *,'WARNING: PROBLEM IN TRIDIAGONAL SYSTEM' PRINT *,'WARNING: ',NMHERE,' TO THE LEFT OF ',I PRINT *,'WARNING: THINGS MAY DETERIORATE' PRINT *,'FIXED: ADJUSTED INTERVALS' ENDIF ENDIF 20 CONTINUE * FIND EIGENVALUES DO 510 I = 1, J IF (DEBUG.GT.1) THEN PRINT *,'UNDONE ',I,UNDONE(I) ENDIF IF (UNDONE(I).EQ.0) THEN IF (I.EQ.1) THEN LEFT = THETA(1) - BJ(1) RIGHT = THETA(1) NTHETA(I) = (LEFT+RIGHT)/2.0D0 NBJ(I) = BJ(I) CALL NEWTON(J,ALPHA,BET2,NTHETA(I),NBJ(I),EPS, C LEFT,RIGHT,DEBUG,I,HUGE) ELSE IF (I.EQ.J) THEN LEFT = THETA(I-1) RIGHT = THETA(I-1) + BJ(I-1) NTHETA(I) = (LEFT+RIGHT)/2.0D0 NBJ(I) = BJ(I) CALL NEWTON(J,ALPHA,BET2,NTHETA(I),NBJ(I),EPS, C LEFT,RIGHT,DEBUG,I,HUGE) ELSE LEFT = THETA(I-1) RIGHT = THETA(I) NTHETA(I) = (THETA(I)+THETA(I-1))/2.0D0 NBJ(I) = BJ(I) CALL NEWTON(J,ALPHA,BET2,NTHETA(I),NBJ(I),EPS, C LEFT,RIGHT,DEBUG,I,HUGE) ENDIF ELSE IF (UNDONE(I).NE.I) THEN ERRTOL = ABS(ANORM*SQRT(EPS)*THETA(I-1)) IF (BJ(I-1).LT.ERRTOL) THEN NBJ(I) = ERRTOL NTHETA(I) = THETA(I-1) ELSE NBJ(I) = BJ(I-1) THETA(I) = THETA(I-1) + (BJ(I-1)/2.0D0) ENDIF LEFT = THETA(I-1) - ABS(THETA(I-1)*EPS*128.0D0) RIGHT = THETA(I-1) + BJ(I-1) ELSE ERRTOL = ABS(ANORM*SQRT(EPS)*THETA(I)) IF (BJ(I).LT.ERRTOL) THEN NBJ(I) = ERRTOL LEFT = THETA(I) - NBJ(I) RIGHT = THETA(I) + ABS(EPS*THETA(I)*128.0D0) ELSE NBJ(I) = BJ(I) LEFT = THETA(I) - BJ(I) RIGHT = THETA(I) + ABS(EPS*THETA(I)*128.0D0) THETA(I) = THETA(I) - (BJ(I)/2.0D0) ENDIF ENDIF CALL NEWTON(J,ALPHA,BET2,NTHETA(I),NBJ(I),EPS, C LEFT,RIGHT,DEBUG,I,HUGE) ENDIF NBJ(I) = NBJ(I)*BETP1 IF (DEBUG.GT.1) THEN PRINT *,'POST = ',NTHETA(I),NBJ(I) ENDIF 510 CONTINUE DO 500 I = 1, J BJ(I) = NBJ(I) THETA(I) = NTHETA(I) 500 CONTINUE ENDIF RETURN END C$FORTRAN EIGMAT *********************************************************************** * LANZ SOFTWARE PACKAGE * * FILENAME: EIGMAT.MSC * * AUTHOR: MARK JONES * * PURPOSE: MATRIX ROUTINES USED IN EIGENVALUE CALCULATIONS * *********************************************************************** * SEVERAL ROUTINES ASSOCIATED WITH EIGENVECTOR * CALCULATIONS * A ROUTINE TO MULTIPLY A TRIDIAGONAL MATRIX, T, BY A VECTOR SUBROUTINE TRIMLT(N,ALPHA,BETA,RESULT,INVEC) * THE DIMENSION OF T INTEGER N * THE DIAGONAL AND SUBDIAGONAL OF T DOUBLE PRECISION ALPHA(*), BETA(*) * RESULT VECTOR DOUBLE PRECISION RESULT(*) * IN VECTOR DOUBLE PRECISION INVEC(*) * INTERNAL VARIABLES * ALL COUNT VARIABLES INTEGER I IF (N.EQ.0) RETURN IF (N.EQ.1) THEN RESULT(1) = ALPHA(1)*INVEC(1) RETURN ENDIF RESULT(1) = ALPHA(1)*INVEC(1) + BETA(2)*INVEC(2) DO 10 I = 2, N-1 RESULT(I) = BETA(I)*INVEC(I-1) + ALPHA(I)*INVEC(I) + C BETA(I+1)*INVEC(I+1) 10 CONTINUE RESULT(N) = BETA(N)*INVEC(N-1) + ALPHA(N)*INVEC(N) RETURN END * A ROUTINE TO MULTIPLY A FULL MATRIX BY A VECTOR SUBROUTINE FLLMLT(LDI,N,M,RESULT,MAT,INVEC,IOPTN) * THE LEADING INDEX INTEGER LDI * THE NUMBER OF ROWS INTEGER N * THE NUMBER OF COLUMNS INTEGER M * THE RESULT VECTOR DOUBLE PRECISION RESULT(*) * THE MATRIX DOUBLE PRECISION MAT(LDI,*) * THE IN VECTOR DOUBLE PRECISION INVEC(*) * MULTIPLY BY MAT OR TRANSPOSE OF MAT INTEGER IOPTN * INTERNAL VARIABLES * ALL COUNT VARIABLES INTEGER I,J IF (IOPTN.EQ.0) THEN DO 5 I = 1, N RESULT(I) = 0.0D0 5 CONTINUE DO 10 I = 1, M DO 20 J = 1, N RESULT(J) = RESULT(J) + MAT(J,I)*INVEC(I) 20 CONTINUE 10 CONTINUE ELSE DO 30 I = 1, M RESULT(I) = 0.0D0 DO 40 J = 1, N RESULT(I) = RESULT(I) + MAT(J,I)*INVEC(J) 40 CONTINUE 30 CONTINUE ENDIF RETURN END * A ROUTINE FOR COMPUTING EIGENVECTORS FROM S AND Q SUBROUTINE CMPVEC(N,Q,QLAST,S,ALPHA,BETA,THETA,Y,PROB) * THE ORDER OF THE MATRIX INTEGER N * THE Q(J) ARRAY DOUBLE PRECISION Q(N,1) * THE LAST INDEX OF THE Q ARRAY INTEGER QLAST * THE S VECTOR DOUBLE PRECISION S(*) * DIAGONAL OF T(J) DOUBLE PRECISION ALPHA(*) * OFF DIAGONAL OF T(J) DOUBLE PRECISION BETA(*) * THE EIGENVALUE DOUBLE PRECISION THETA * THE EIGENVECTOR DOUBLE PRECISION Y(*) * THE TYPE OF TRANSFORMATION INTEGER PROB * THE FOLLOWING ARE INTERNAL VARIABLES * COMPUTE THE ADJUSTED VECTOR, W, FROM S AS SUGGESTED BY PARLETT * W(1) = ALPHA(1)*S(1) + BETA(2)*S(2) * DO 5 I = 2, QLAST-1 * W(I) = BETA(I)*S(I-1) + ALPHA(I)*S(I) + * C BETA(I+1)*S(I+1) *5 CONTINUE * IF (QLAST.GT.1) THEN * W(QLAST) = BETA(QLAST)*S(QLAST-1) + * C ALPHA(QLAST)*S(QLAST) * ENDIF * W(QLAST+1) = S(QLAST) * BETA(QLAST+1) * DO 8 I = 1, QLAST+1 * W(I) = W(I) / THETA *8 CONTINUE IF (PROB.EQ.1) THEN S(QLAST+1) = (S(QLAST) * BETA(QLAST+1))/(THETA-1) ELSE S(QLAST+1) = (S(QLAST) * BETA(QLAST+1))/THETA ENDIF * DO 50 I = 1, QLAST+1 * PRINT *,'S = ',S(I) *50 CONTINUE * NOW COMPUTE Y FROM Q*W CALL FLLMLT(N,N,QLAST+1,Y,Q,S,0) RETURN END C$FORTRAN NWSHFT *********************************************************************** * LANZ SOFTWARE PACKAGE * * FILENAME: FNDSHT.MSC * * AUTHOR: MARK JONES * * PURPOSE: FINDS A NEW SHIFT(S) FOR LANZ * *********************************************************************** * PREVIOUS TO CALLING THE ROUTINE, WE MAY HAVE A LEFT SHIFT AND * A RIGHT SHIFT AND THE PURPOSE OF THIS ROUTINE IS TO FIND GOOD * NEW LEFT AND RIGHT SHIFTS. IT IS POSSIBLE TO NOT HAVE A NEW * SHIFT AT EITHER OR BOTH POSITIONS. THE INFORMATION USED IF * FROM THE OLD EIGENVALUES AND THE LATEST RUNS (R1 AND R2) * AT THE LEFT AND RIGHT SHIFTS SUBROUTINE NWSHFT(J,THETA,BJ,R1J,R1THET,R1BJ,R2J,R2THET,R2BJ, C LFTVAL,RGHTVL,NWLEFT,NWRGHT,DEBUG,KNORM,MNORM, C EPS,DELTA,OSIGMA) * THE ORDER OF THE TRIDIAGONAL MATRICES INTEGER J, R1J, R2J * THE EIGENVALUES FOUND SO FAR DOUBLE PRECISION THETA(*), R1THET(*), R2THET(*) * THE ERROR BOUNDS ON THESE EIGENVALUES DOUBLE PRECISION BJ(*), R1BJ(*), R2BJ(*) * CURRENT LEFT EDGE AND CURRENT RIGHT EDGE DOUBLE PRECISION LFTVAL, RGHTVL * DID WE UPDATE THE LEFT OR RIGHT EDGES LOGICAL NWLEFT, NWRGHT * THE DEBUGGING LEVEL INTEGER DEBUG * THE NORMS OF K AND M DOUBLE PRECISION KNORM, MNORM * THE MACHINE EPS DOUBLE PRECISION EPS * THE DESIRED ACCURACY FOR THE EIGENVALUES DOUBLE PRECISION DELTA * THE ORIGINAL SIGMA DOUBLE PRECISION OSIGMA * SOME INTERNAL VARIABLES * CONVERGED BOUNDARIES DOUBLE PRECISION CLEFT, CRIGHT * UNCONVERGED BOUNDARIES DOUBLE PRECISION UNCLFT, UNCRGT * ERROR BOUNDS FOR ABOVE DOUBLE PRECISION UNCLBJ, UNCRBJ * FOUND VALUES TO THE LEFT OR RIGHT LOGICAL FLEFT, FRIGHT * THE MINIMUM DISTANCE THAT A SHIFT MAY BE FROM A POSSIBLE * EIGENVALUE DOUBLE PRECISION MNDIST * THE MAXIMUM PERCENTAGE OF THE EIGENVALUE THAT THE ERROR BOUND * MAY BE, BEFORE WE IGNORE IT DOUBLE PRECISION FRACT * THE OLD VALUES OF PARAMETERS DESCRIBED ABOVE LOGICAL ONWRGT, ONWLFT DOUBLE PRECISION ORGTVL, OLFTVL * SOME PARAMETERS USED TO DETERMINE WHERE TO TAKE A SHIFT DOUBLE PRECISION DELERR, ENU, ENORM EXTERNAL CHKBND, ADJSHT * FIND WHERE UNCONVERGED EIGENVALUE BOUNDARIES ARE FRACT = 0.10D0 CLEFT = LFTVAL CRIGHT = RGHTVL UNCLFT = LFTVAL UNCRGT = RGHTVL UNCLBJ = 0.0D0 UNCRBJ = 0.0D0 FLEFT = .FALSE. FRIGHT = .FALSE. CALL CHKBND(J,THETA,BJ,CLEFT,CRIGHT,UNCLFT,UNCRGT, C UNCLBJ,UNCRBJ,LFTVAL,RGHTVL,FRACT,OSIGMA,FLE