** Sample main program ... remove the asterisks in the first column ** to obtain a working code ** ********************************************************************* ** * PROGRAM MAIN ** ********************************************************************* ** * PARAMETER (NW=150000,NH=10000) * INTEGER IWRK(NW) * REAL HWRK(NH) **------------------------------------------------------------------- ** DEC20 CODE **------------------------------------------------------------------- ** write (*,*) ' give ibase ' ** read(*,*) ibase ** IBASE = 100000 ** LGINT = 34359738367 ** ITRMI = 5 ** ITRMO = 6 **------------------------------------------------------------------- ** CRAY CODE **------------------------------------------------------------------- ** CALL LINK(' UNIT5=TERMINAL, UNIT10 = D, UNIT11 = R11, ** * UNIT12 = (R12, CREATE, HARDCOPY), ** * UNIT13 = (R13, CREATE, HARDCOPY) //') ** LGINT = 9223372036854775807 ** IBASE = 10000000 ** ITRMI = 5 ** ITRMO = 5 **------------------------------------------------------------------- ** VAX CODE **------------------------------------------------------------------- * LGINT = 2147483647 * IBASE = 10000 * ITRMI = 5 * ITRMO = 6 **------------------------------------------------------------------- ** SUN CODE **------------------------------------------------------------------- ** LGINT = 2147483647 ** IBASE = 10000 ** ITRMI = 5 ** ITRMO = 6 **------------------------------------------------------------------- ** MACINTOSH CODE **------------------------------------------------------------------- ** LGINT = 2147483647 ** IBASE = 10000 ** ITRMI = 9 ** ITRMO = 9 **------------------------------------------------------------------- * CALL SECOND(T1) * CALL DIALOG(IWRK,NW,HWRK,NH,ITRMI,ITRMO,IBASE,LGINT) * CALL SECOND(T2) * WRITE(ITRMO,*) ' TOTAL TIME = ',T2-T1 * STOP * END SUBROUTINE GOLIA(IOUT) * THIS SUBROUTINE WRITES AN IDENTIFYING MESSAGE TO CHANNEL IOUT. * THE MESSAGE SHOULD BE MODIFIED APPROPRIATELY WHENEVER THE PROGRAM * IS MODIFIED. IF NO MESSAGE IS DESIRED THE FOLLOWING WRITE STATEMENT * SHOULD BE COMMENTED OUT. WRITE (IOUT,100) 100 FORMAT(' >>> this is GOLIATH, version 2.03, April 12, 1989') RETURN END * * * If the following dummy routine is replaced by a subroutine that * assignes to T the amount of CPU time elapsed since the program started * then the output will contain some timing data. * SUBROUTINE SECOND(T) REAL T RETURN END * SUBROUTINE ANALYZ(ITRMO,IBASE,LGINT,IORAT,IOINT,IORES,IOLOG,LOG, X JOB,ITYPE,ISTOP,IBLANK,IWRK,NWRK,HWRK,NH) *----------------------------------------------------------------------- * *** Goliath *** * * A system for the exact analysis of sparse rectangular rational * linear systems. Developed by * * Peter Alfeld and David Eyre * Department of Mathematics * University of Utah * Salt Lake City, Utah 84112 * (801) 581-6842 or (801) 581-6851 * ALFELD@SCIENCE.UTAH.EDU * EYRE@SCIENCE.UTAH.EDU * * This system obtains a linear system with rational or integer coefficients * from a file, and prints the results of the analysis into a file. * * For a detailed documentation see the Collected Algorithms of the ACM * * Most of the input parameters can be set interactively using the * SUBROUTINE DIALOG. The parameters have the following meanings: * * ITRMO: output channel number, normally the terminal. * * IBASE: The base of the fixed radix output, normally a power of 10. * It must be possible to represent the square of IBASE in one * INTEGER word. * * LGINT: The largest integer that can be represented on the machine * * IORAT: Channel number of a rational input file * * IOINT: Channel number of an integer intermediate or input file. * * IORES: Channel number for the result file * * IOLOG: Channel number for file receiving intermediate output. * * LOG: Amount of intermediate output, directed to IOLOG. * Possibilities include: * 0 no output * 1 Operation statistics at end of run * 2 Elimination info during first run only * 3 All matrix information, pass 1 only * 4 All info, all passes, caution! * * JOB: Type of task to be completed. Possibilities include: * 0 Analysis of Rank and Row dependencies * 1 Solution of Linear Systems * 2 Analysis of Null Space * 3 Solution of Linear Systems and Analysis of Null Space * 4 Conversion of rational file to integer file * * ITYPE: Type of input file. Possibilities include: * 0 Arbitrary Precision Integer * 1 Single Precision Rational * 2 Arbitrary Precision Rational * * ISTOP: stopping criterion. Possibilities include: * -2 Modified Hadamard bound * -1 Convergence criterion * 0 Hadamard bound * N>0 specified number of passes * * IBLANK: describes form of fixed radix output * 0 no blanks between terms * 1 blank between spaces * * IWRK: 1 dimensional Integer work array of sufficient size * * NWRK: Size of IWRK * * HWRK: 1 dimensional real work array. * * NH: size of HWRK. It must be greater than * # of rows + # of columns + # of right hand sides * * INTEGER IPARM(0:15),IWRK(NWRK),MPARM(7),MINFO(4) REAL HWRK(NH) *----------------------------------------------------------------------- * Transfer parameters to IPARM and MINFO * These arrays are explained in the subroutine RSDANA *----------------------------------------------------------------------- IF (ITYPE.GT.0) THEN READ(IORAT,*) (MINFO(I),I=1,4) REWIND IORAT ELSEIF (ITYPE.EQ.0) THEN READ(IOINT,*) (MINFO(I),I=1,4) REWIND IOINT ENDIF IPARM(0) = MINFO(1) IPARM(1) = MINFO(2) + MINFO(3) IPARM(2) = MINFO(3) IPARM(3) = IOINT IPARM(4) = IORES IPARM(5) = IOLOG IPARM(6) = ITRMO IPARM(7) = LOG IPARM(8) = JOB IPARM(9) = IBASE IPARM(10) = ISTOP IPARM(11) = LGINT IPARM(12) = 0 IPARM(13) = IBLANK IPARM(14) = ITYPE IPARM(15) = IORAT *------------------------------------------------------------------- * EXECUTE MPCON *------------------------------------------------------------------- MPARM(7) = 0 IF (IPARM(14) .GT. 0) THEN MPARM(1) = IPARM(15) MPARM(2) = IPARM(3) MPARM(3) = IPARM(9) MPARM(4) = 0 IF ( IPARM(14) .EQ. 1 ) MPARM(4) = 1 MPARM(5) = IPARM(7) MPARM(6) = IPARM(6) MPARM(7) = INT(SQRT(FLOAT(IPARM(11)))/1.3) CALL MPCON (MPARM,IWRK,NWRK) ENDIF IF (IPARM(8).NE.4) THEN IF ( IPARM(10) .GT. 0 ) THEN IPARM(12) = IPARM(12) + 100 ELSEIF ( MPARM(7) .GT. 0 ) THEN IPARM(12) = MPARM(7) + 10 ELSE J = (NWRK-1) / 2 IPARM(12) = IPR(IPARM,IWRK(1),IWRK(J+1),J) + 10 ENDIF *------------------------------------------------------------------- * EXECUTE RSDANA *------------------------------------------------------------------- JWRK = IWORK (IPARM) NSTOR = MAX0(0,(NWRK-JWRK-2)/2) CALL RSDANA (IPARM,IWRK(JWRK+1),NSTOR,IWRK(1),JWRK,HWRK,NH) ENDIF RETURN END * * SUBROUTINE DIALOG(IWRK,NWRK,HWRK,NH,ITRMI,ITRMO,IBASE,LGINT) * * THIS SUBROUTINE ACCESSES ANALYZ AFTER A DIALOG WITH THE USER. * DATA ARE READ FROM A FILE, AND RESULTS ARE DIRECTED TO A FILE. * * THE PARAMETERS OF DIALOG HAVE THE FOLLOWING MEANINGS: * * IWRK: A 1-DIMENSIONAL INTEGER WORK ARRAY. THE LARGER IT IS THE LAR * ARE THE PROBLEMS THAT CAN BE SOLVED. HOWEVER, IT SHOULD NOT * EXCEED THE MAXIMUM REQUIRED WORKSPACE BY TOO LARGE A FACTOR * THEN DISK OPERATIONS REQUIRED BY SWAPPING MAY DECREASE EFFICI * NWRK: PHYSICAL DIMENSION OF IWRK. * HWRK: REAL WORK ARRAY * NH: SIZE OF HWRK, REQUIRED MINIMUM VALUE IS * NH > (# OF ROWS) + (# OF COLUMNS) + (# OF RHS) * ITRMI: UNIT NUMBER FOR INPUT FROM TERMINAL * ITRMO: UNIT NUMBER FOR OUTPUT TO TERMINAL * IBASE: THE BASE FOR MULTIPLE PRECISION OUTPUT. FOR BETTER LEGIBILTI * IT SHOULD BE A POWER OF 10. * LGINT: THE LARGEST INTEGER THAT CAN BE EXPRESSED ON THE MACHINE. * * * * THE PROBLEM SOLVED IS: * * AX = B * WHERE: * A IS AN NXM RATIONAL MATRIX * X IS THE MXP RATIONAL UNKNOWN MATRIX * B IS AN NXP RATIONAL RIGHT HAND SIDE. * * INPUT FOR THIS PROGRAM IS INTERACTIVE, AND IT WILL ACCEPT * FILES IN EITHER RATIONAL OR INTEGER FORMAT. RATIONAL FORMAT * FILES ARE PROCESSED BY MPCON. THEN THE ACTUAL ANALYSIS IS * PERFORMED IN RSDANA. * * THIS SUBROUTINE USES IO CHANNELS 90,91,92,93 INTERNALLY * AS FOLLOWS: * I90: INTEGER INPUT FILE * I91: INTEGER FILE NAME * I92: RESULT FILE NAME * I93: LOG FILE NAME * IF OTHER CHANNEL NUMBERS ARE REQUIRED THE FOLLOWING PARAMETER * STATEMENT SHOULD BE ALTERED APPROPRIATELY: * PARAMETER(I90=90,I91=91,I92=92,I93=93) * DATA STRUCTURE DECLARATION *------------------------------------------------------------------- INTEGER IWRK(NWRK) REAL HWRK(NH) CHARACTER*25 FNM1,FNM2 CALL GOLIA(ITRMO) *------------------------------------------------------------------- * INPUT FROM THE USER - DATA FILE CONTAINING MATRIX AND RHS *------------------------------------------------------------------- WRITE (ITRMO,1100) NWRK READ (ITRMI,*) NW IF (0.GE.NW.OR.NW.GT.NWRK) NW = NWRK WRITE (ITRMO,1080) READ (ITRMI,*) ITYPE IF ( ITYPE .GT. 0 ) THEN WRITE (ITRMO,1010) READ (ITRMI,1000) FNM1 OPEN (I91,FILE=FNM1,STATUS='OLD') WRITE (ITRMO,1090) READ (ITRMI,1000) FNM1 OPEN (I90,FILE=FNM1,STATUS='UNKNOWN') ELSEIF (ITYPE .EQ. 0 ) THEN WRITE (ITRMO,1010) READ (ITRMI,1000) FNM1 OPEN (I90,FILE=FNM1,STATUS='OLD') ENDIF IOINT = I90 IORAT = I91 *------------------------------------------------------------------- * INPUT FROM THE USER - DATA OUTPUT *------------------------------------------------------------------- WRITE (ITRMO,1030) READ (ITRMI,*) JOB IF (JOB.LT.4) THEN *------------------------------------------------------------------- * INPUT FROM THE USER - OUTPUT FILE THAT WILL STORE RESULTS *------------------------------------------------------------------- WRITE (ITRMO,1020) READ (ITRMI,1000) FNM2 IF (FNM2(1:1) .EQ. '0') THEN IORES = 0 ELSE OPEN (I92,FILE=FNM2,STATUS='UNKNOWN') IORES = I92 ENDIF WRITE (ITRMO,1040) READ (ITRMI,*) IBLANK *------------------------------------------------------------------- * INPUT FROM THE USER - TERMINATION CRITERION *------------------------------------------------------------------- WRITE (ITRMO,1050) READ (ITRMI,*) ISTOP *------------------------------------------------------------------- * INPUT FROM THE USER - DATA FILE WITH INTERMEDIATE OUTPUT *------------------------------------------------------------------- WRITE (ITRMO,1060) READ (ITRMI,*) LOG IF ( LOG .GT. 1 ) THEN WRITE (ITRMO,1070) READ (ITRMI,1000) FNM2 OPEN (I93,FILE=FNM2,STATUS='UNKNOWN') ENDIF IOLOG = I93 ELSE LOG = 0 ENDIF CALL ANALYZ(ITRMO,IBASE,LGINT,IORAT,IOINT,IORES,IOLOG,LOG, X JOB,ITYPE,ISTOP,IBLANK,IWRK,NW,HWRK,NH) *------------------------------------------------------------------- * CLOSE OPENED FILES AND SHUT DOWN *------------------------------------------------------------------- CLOSE (I90) CLOSE (I91) CLOSE (I92) CLOSE (I93) RETURN 1000 FORMAT (A25) 1010 FORMAT (/' enter input file name: ') 1020 FORMAT (/' enter output file name (enter 0 for no hardcopy):') 1030 FORMAT (/' enter desired computation:',/ * ' 0 = rank, dependencies',/ * ' 1 = rank and rhs solutions',/ * ' 2 = rank and null space basis',/ * ' 3 = rank, rhs and null space basis',/ * ' 4 = rational ----> integer conversion') 1040 FORMAT (/' enter the mp output style:',/ * ' 0 = no spaces between terms',/ * ' 1 = one space between terms') 1050 FORMAT (/' enter termination criterion:',/ * ' -2 = modified hadamard bound',/ * ' -1 = convergence criterion',/ * ' 0 = hadamard bound',/ * ' n = specific pass count') 1060 FORMAT (/' enter intermediate output parameter:',/ * ' 0 = none ',/ * ' 1 = terminal progress report',/ * ' 2 = reduction stage info, pass one only',/ * ' 3 = all matrix, pass one only',/ * ' 4 = all info - all passes ***caution***') 1070 FORMAT (/' enter intermediate output file name: ') 1080 FORMAT (/' enter input file type: ',/ * ' 0 = integer ',/ * ' 1 = rational',/ * ' 2 = arbitrary precision rational') 1090 FORMAT (/' enter integer file name: ') 1100 FORMAT (/' enter size of workspace (up to',I12,')') END * ******************************************************************** * SUBROUTINE GJPIV (K,JZRO,IWRK,NWRK,IHTBL,NSTOR) * * PERFORMS COLUMN PIVOTING OF THE AUGMENTED MATRIX. THE PIVOT * SELECTION IS BASED ON * * 1. AKJ NOT EQUAL TO ZERO * 2. JTH COLUMN HAVING THE LEAST NON ZERO ELEMENTS * * ROW PIVOTING IS PERFORMED IF THE KTH ROW IS THE 0 * VECTOR BY ROTATING THE KTH ROW THE LAST ROW IN THE * MATRIX. THEN THE K+1ST ROW BECOMES THE KTH ROW AND * COLUMN PIVOTING IS PERFORMED ON THE IT AS ABOVE. * THE NUMBER OF NON-ZERO ELEMENTS IN ANY GIVEN COLUMN * IS COUNTED WITH THE IPIV ARRAY. * * GLOBAL VARIABLES: * * 1. K - INTEGER - CURRENT STAGE * 2. JZRO - INTEGER - SIGNALS WHEN MATRIX HAS ONLY ZEROS * BELOW THE KTH ROW, I.E. REDUCTION IS COMPLETE * = 0 THEN ONLY ZEROS LEFT * = 1 THEN CONTINUE * 3. IWRK - NWRK INTEGER ARRAY - SEE DRIVER FOR SPECIFICATION * 4. IHTBL - NSTOR INTEGER ARRAY - SEE DRIVER FOR SPECIFICATION * * COMMON VARIABLES CHANGED: * * LOCAL VARIABLES: * * 1. NKROW - INTEGER - FLAG SIGNALING ROW PIVOTING NECESSARY * 2. IPCOL - INTEGER - POTENTIAL PIVOT COLUMN * 3. IDUM - INTEGER - VOLATILE * * NOTE: COMMON VARIABLES LISTED IN DRIVER * ******************************************************************** * INTEGER IHTBL(2,0:NSTOR),IWRK(NWRK) INTEGER IDUM,IPCOL,NKROW COMMON /WPARM/ MINFO,IHPRM,IROW,ICOL,IDEP,IPIV,LST,IPRM,IPIN, * IAPRM,KODE,IRHS,JWRK,KNUM,KGCD,IMAT,KUN13,IBASE, * ITER,ICPRM,IDIAG,ITRM,IBASIC,ITRI,INULL,IRANK,IZERO,LMAX, * NPRME,MHZ,ISOLH,KUN12 *------------------------------------------------------------------- * INITIALIZE - CHECK TO SEE IF MATRIX IS REDUCED *------------------------------------------------------------------- 10 CONTINUE NKROW = 0 IMIN = IWRK(MINFO+1) + 1 IF ( IWRK(IROW+K) .EQ. IWRK(IDEP+2) ) THEN JZRO = 0 IMIN = 0 GO TO 40 ENDIF *------------------------------------------------------------------- * FIND NONZERO COLUMN ENTRIES IN ROW K. IF A NONZERO PIVOT * IS FOUND SEARCH ITS COLUMN. ELSE CHECK TO SEE IF IT IS A * NEW MINIMUM AND CONTINUE. *------------------------------------------------------------------- DO 20, J = K, IWRK(MINFO+2) KEY = IRHASH(IWRK(IROW+K),IWRK(ICOL+J),0,IWRK(IHPRM+1), * IHTBL,NSTOR) IF ( IHTBL(2,KEY) .NE. 0 ) THEN NKROW = NKROW + 1 IDUM = IWRK(IPIV+IWRK(ICOL+J)) IF (IDUM .LT. IMIN) THEN IMIN = IDUM IPCOL = J ENDIF ENDIF 20 CONTINUE *------------------------------------------------------------------- * IF NO ENTRY IN THE K ROW IS NONZERO THEN ROW PIVOT AND * REPEAT THE COLUMN MINIMIZATION. *------------------------------------------------------------------- IF (NKROW .EQ. 0) THEN IWRK(IDEP+1) = IWRK(IDEP+1) + 1 IWRK(IDEP+IWRK(IDEP+1)+1) = IWRK(IROW+K) IDUM = IWRK(IROW+K) DO 30, I = K+1, IWRK(MINFO+1) IWRK(IROW+I-1) = IWRK(IROW+I) 30 CONTINUE IWRK(IROW+IWRK(MINFO+1)) = IDUM GO TO 10 ENDIF *------------------------------------------------------------------- * SWAP ROWS FROM MINIMIZING COUNTER *------------------------------------------------------------------- IDUM = IWRK(ICOL+K) IWRK(ICOL+K) = IWRK(ICOL+IPCOL) IWRK(ICOL+IPCOL) = IDUM 40 CONTINUE RETURN END * ******************************************************************** * SUBROUTINE GJRES (IWRK,NWRK,IHTBL,NSTOR) * * COMPUTES INTERMEDIATE RESULTS MODULO ICPRM OF THE PRINCIPAL * MINOR, PARTICULAR SOLUTION, AND NULL SPACE BY RECOVERING, * FOR EACH ELEMENT IN THE SOLUTION, THE NEW RESIDUES AND ITS * PREVIOUS MIXED RADIX FORM AND UPDATES THE MIXED RADIX FORM * IN RADIXM. CHECKS THE RANK OF THE CURRENT PASS AND CHECKS * THE RIGHT HAND SIDE FOR CONSISTENCY. HASHES SOLUTION INTO * SOLUTION HASH TABLE. * * GLOBAL VARIABLES: * * 1. IWRK - NWRK INTEGER ARRAY - SEE DRIVER FOR SPECIFICATION * 2. IHTBL - NSTOR INTEGER ARRAY - SEE DRIVER FOR SPECIFICATION * * COMMON VARIABLES USED: * 1. IRANK - RANK FLAG THAT IS SET * = 1 RANK INCREASED * = 0 RANK UNCHANGED * = -1 RANK DECREASED * 2. IZERO - OUTPUT FLAG FOR ZERO COEFFIENTS * = 0 IF ALL COEFFICIENTS ARE ZERO * = 1 IF ALL COEFFICIENTS ARE NON-ZERO * * LOCAL VARIABLES: * * 1. JZ - INTEGER - VOLATILE * ******************************************************************** * INTEGER IHTBL(2,0:NSTOR),IWRK(NWRK) INTEGER JZ LOGICAL TEST COMMON /WPARM/ MINFO,IHPRM,IROW,ICOL,IDEP,IPIV,LST,IPRM,IPIN, * IAPRM,KODE,IRHS,JWRK,KNUM,KGCD,IMAT,KUN13,IBASE, * ITER,ICPRM,IDIAG,ITRM,IBASIC,ITRI,INULL,IRANK,IZERO,LMAX, * NPRME,MHZ,ISOLH,KUN12 *------------------------------------------------------------------- * FROM THE MAIN DIAGONAL CHECK FOR DECREASED RANK. *------------------------------------------------------------------- IZERO = 1 IRANK = 0 JZ = 1 IF ( ITER .NE. 1 ) THEN DO 10, I = 1, IWRK(MINFO+7) KEY = IRHASH(IWRK(IROW+I),IWRK(ICOL+I),0,IWRK(IHPRM+1), * IHTBL,NSTOR) JZ = JZ * IHTBL(2,KEY) 10 CONTINUE IF ( JZ .NE. 1 ) THEN IRANK = -1 RETURN ENDIF *------------------------------------------------------------------- * CHECK FOR INCREASED RANK AND SOLUTION CONSISTENCY WITH THE ROW * COUNTER VECTOR *------------------------------------------------------------------- DO 30, I = IWRK(MINFO+7)+1, IWRK(MINFO+1) JZ = IWRK(IROW+I) IF ( IWRK(IPIV+JZ) .NE. 0 ) THEN DO 20, J = 1, IWRK(MINFO+3) KEY = IRHASH(JZ,IWRK(MINFO+2)+J,0,IWRK(IHPRM+1), * IHTBL,NSTOR) IF ( IHTBL(2,KEY) .NE. 0 ) THEN IWRK(IRHS+J) = 1 IWRK(IPIV+JZ) = IWRK(IPIV+JZ) - 1 ENDIF 20 CONTINUE IF ( IWRK(IPIV+JZ) .NE. 0 ) IRANK = 1 ENDIF 30 CONTINUE ENDIF *------------------------------------------------------------------- * UPDATE COEFFICIENT OF THE LEADING MINOR *------------------------------------------------------------------- TEST = .TRUE. CALL RADIXM (KNUM,IDIAG,TEST,IWRK,NWRK) *------------------------------------------------------------------- * COMPUTE THE PARTICULAR SOLUTION FOR THE CURRENT PASS WHEN * GAUSS-JORDAN DECOMPOSTION IS USED. TAKE INTO ACCOUNT * PIVOTING. *------------------------------------------------------------------- ITN = ITER-1 IF (IBASIC.EQ.1 .AND. ITRI.EQ.0) THEN DO 50, J = IWRK(MINFO+2)+1, IWRK(MINFO+3)+IWRK(MINFO+2) JS = J DO 40, I = 1, IWRK(MINFO+7) IST = IWRK(IROW+I) IS = IWRK(ICOL+I) KEY = IRHASH(IST,JS,0,IWRK(IHPRM+1),IHTBL,NSTOR) KK = MOD(IHTBL(2,KEY)*IDIAG,ICPRM) CALL JRHASH (IS,JS,LST,ITN,0,JZ,IWRK,NWRK,IHTBL,NSTOR) IF ( KK.NE.0 .OR. JZ.NE.0 ) THEN CALL RADIXM (LST,KK,TEST,IWRK,NWRK) CALL JSHASH (IS,JS,LST,ITER,IWRK,NWRK,IHTBL,NSTOR) ENDIF 40 CONTINUE 50 CONTINUE *------------------------------------------------------------------- * COMPUTE AND OUTPUT THE PARTICULAR SOLUTION FOR THE CURRENT PASS * WHEN TRIANGULAR DECOMPOSTION IS USED - BACK SUBSTITUTION *------------------------------------------------------------------- ELSEIF (IBASIC.EQ.1 .AND. ITRI.EQ.1) THEN DO 80, J = IWRK(MINFO+2)+1, IWRK(MINFO+3)+IWRK(MINFO+2) JS = J DO 70, I = IWRK(MINFO+7), 1, -1 IST = IWRK(IROW+I) IS = IWRK(ICOL+I) KEY = IRHASH(IST,JS,0,IWRK(IHPRM+1),IHTBL,NSTOR) KK = IHTBL(2,KEY) DO 60, K = I+1, IWRK(MINFO+7), 1 KEY = IRHASH(IST,IWRK(ICOL+K),0,IWRK(IHPRM+1), * IHTBL,NSTOR) KK = MOD( KK - MOD(IWRK(KODE+K)*IHTBL(2,KEY),ICPRM), * ICPRM) 60 CONTINUE IWRK(KODE+I) = KK KK = MOD(KK*IDIAG,ICPRM) CALL JRHASH (IS,JS,LST,ITN,0,JZ,IWRK,NWRK,IHTBL,NSTOR) IF ( KK.NE.0 .OR. JZ.NE.0 ) THEN CALL RADIXM (LST,KK,TEST,IWRK,NWRK) CALL JSHASH (IS,JS,LST,ITER,IWRK,NWRK,IHTBL,NSTOR) ENDIF 70 CONTINUE 80 CONTINUE ENDIF *------------------------------------------------------------------- * OUTPUT NULL SPACE FOR THE CURRENT PASS *------------------------------------------------------------------- IF (INULL .EQ. 1) THEN DO 100, I = 1, IWRK(MINFO+7) IS = IWRK(ICOL+I) DO 90, J = IWRK(MINFO+7)+1, IWRK(MINFO+2) JS = J KEY = IRHASH(IWRK(IROW+I),IWRK(ICOL+J),0,IWRK(IHPRM+1), * IHTBL,NSTOR) KK = MOD(IHTBL(2,KEY)*IDIAG,ICPRM) CALL JRHASH (IS,JS,LST,ITN,0,JZ,IWRK,NWRK,IHTBL,NSTOR) IF ( KK.NE.0 .OR. JZ.NE.0 ) THEN CALL RADIXM (LST,KK,TEST,IWRK,NWRK) CALL JSHASH (IS,JS,LST,ITER,IWRK,NWRK,IHTBL,NSTOR) ENDIF 90 CONTINUE 100 CONTINUE ENDIF IF ( TEST ) IZERO = 0 RETURN END * ******************************************************************** * SUBROUTINE GJRRED (K,IP,IWRK,NWRK,IHTBL,NSTOR) * * ROW REDUCES THE KTH COLUMN OF THE AUGMENTED MATRIX MODULO * ICPRM WITH GAUSS OR GAUSS-JORDAN ROW REDUCTION AS APPROP. * PIVOTING FOR THE KTH STAGE SHOULD ALREADY BE DECIDED. * ASSUMES THAT IF AKK IS ZERO THEN THE MATRIX IS REDUCED TO * IDENTITY FORM AND RETURNS. COUNTS THE NUMBER OF NON-ZERO * ELEMENTS IN ANY GIVEN COLUMN ON THE FIRST PASS OR ROW ON * SUCCESSIVE PASSES IN THE IPIV ARRAY. THIS IS USED IN * PIVOTING ON THE FIRST PASS AND RANK DETERMINATION ON * SUCCESSIVE PASSES. * * GLOBAL VARIABLES: * * 1. K - INTEGER - CURRENT STAGE * 2. IP - INTEGER - FLAG FOR PIVOT COUNTER * = 0 COUNT COLUMN ZEROES FOR PIVOTING * = ANYTHING ELSE THEN COUNT ROW ZEORES FOR * RANK DETERMINATION * 3. IWRK - NWRK INTEGER ARRAY - SEE DRIVER FOR SPECIFICATION * 4. IHTBL - NSTOR INTEGER ARRAY - SEE DRIVER FOR SPECIFICATION * * COMMON VARIABLES CHANGED: * * LOCAL VARIABLES: * * 1. JAIJ,JAKK,JAIK - INTEGERS - OBVIOUS MATRIX ELEMENTS * 2. IAIJ - INTEGER - AIJ BEFORE REDUCTION FOR PIVOT COUNTER * 3. KK - VOLATILE * ******************************************************************** * INTEGER IHTBL(2,0:NSTOR),IWRK(NWRK) INTEGER JAIJ,JAKK,JAIK,IAIJ COMMON /WPARM/ MINFO,IHPRM,IROW,ICOL,IDEP,IPIV,LST,IPRM,IPIN, * IAPRM,KODE,IRHS,JWRK,KNUM,KGCD,IMAT,KUN13,IBASE, * ITER,ICPRM,IDIAG,ITRM,IBASIC,ITRI,INULL,IRANK,IZERO,LMAX, * NPRME,MHZ,ISOLH,KUN12 *------------------------------------------------------------------- * LOAD AKK AND COMPUTE ITS MULTIPICATIVE INVERSE FOR GAUSSIAN * REDUCTION TO THE APPRPRIATE FORM. UPDATE IDIAG WHICH STORES * THE PIVOT ELEMENT INFORMATION. *------------------------------------------------------------------- KEY = IRHASH(IWRK(IROW+K),IWRK(ICOL+K),0,IWRK(IHPRM+1), * IHTBL,NSTOR) JAKK = IHTBL(2,KEY) IF ( JAKK .EQ. 0 ) RETURN IHTBL(2,KEY) = 1 IDIAG = MOD(IDIAG*JAKK,ICPRM) JAKK = INV(JAKK,ICPRM) *------------------------------------------------------------------- * LOAD KTH ROW ARRAY FOR ELIMINATION, MULTIPLY BY JAKK INVERSE * AND RESTORE *------------------------------------------------------------------- DO 10, J = K+1, IWRK(MINFO+2)+IWRK(MINFO+3) KEY = IRHASH(IWRK(IROW+K),IWRK(ICOL+J),0,IWRK(IHPRM+1), * IHTBL,NSTOR) IWRK(LST+J) = MOD(IHTBL(2,KEY)*JAKK,ICPRM) IHTBL(2,KEY) = IWRK(LST+J) 10 CONTINUE *------------------------------------------------------------------- * SEARCH KTH COLUMN FOR NONZERO ENTRIES - PERFORM GAUSSIAN ELIM. * IF NULL SPACE IS NOT NEEDED, ELSE PERFORM GAUSS-JORDAN * ELIMINATION SO THAT * AIJ := AIJ - AIK*AKJ. * NOTE THAT AKJ IS STORED IN IWRK(LST+J). *------------------------------------------------------------------- IF ( ITRI .EQ. 1 ) THEN II = K + 1 ELSE II = 1 ENDIF DO 30, I = II, IWRK(MINFO+1) IF ( I .EQ. K ) GO TO 30 KEY = IRHASH(IWRK(IROW+I),IWRK(ICOL+K),0,IWRK(IHPRM+1), * IHTBL,NSTOR) JAIK = IHTBL(2,KEY) IF (JAIK .NE. 0) THEN IHTBL(2,KEY) = 0 IF ( IP .NE. 0 ) * IWRK(IPIV+IWRK(IROW+I)) = IWRK(IPIV+IWRK(IROW+I)) - 1 DO 20, J = K+1, IWRK(MINFO+2)+IWRK(MINFO+3) KEY = IRHASH(IWRK(IROW+I),IWRK(ICOL+J),0,IWRK(IHPRM+1), * IHTBL,NSTOR) JAIJ = IHTBL(2,KEY) IAIJ = JAIJ IF (JAIJ .EQ. 0) KEY = ISHASH(IWRK(IROW+I),IWRK(ICOL+J), * 0,IWRK(IHPRM+1),IHTBL,NSTOR) JAIJ = MOD( JAIJ - MOD(JAIK*IWRK(LST+J),ICPRM), ICPRM) IHTBL(2,KEY) = JAIJ *------------------------------------------------------------------- * RECORD FILLIN FOR PIVOTING ON FIRST PASS *------------------------------------------------------------------- IF ( IP .EQ. 0 ) THEN IF ( IAIJ .EQ. 0 .AND. JAIJ .NE. 0) THEN IWRK(IPIV+IWRK(ICOL+J)) = IWRK(IPIV+IWRK(ICOL+J)) + 1 IF ( J .LE. IWRK(MINFO+2) ) * IWRK(MINFO+9) = IWRK(MINFO+9) + 1 ELSEIF ( IAIJ .NE. 0 .AND. JAIJ .EQ. 0) THEN IWRK(IPIV+IWRK(ICOL+J)) = IWRK(IPIV+IWRK(ICOL+J)) - 1 ENDIF *------------------------------------------------------------------- * RECORD FILLIN FOR RANK DETERMINATION ON SUBSEQUENT PASSES *------------------------------------------------------------------- ELSE IF ( IAIJ .EQ. 0 .AND. JAIJ .NE. 0) THEN IWRK(IPIV+IWRK(IROW+I)) = IWRK(IPIV+IWRK(IROW+I)) + 1 ELSEIF ( IAIJ .NE. 0 .AND. JAIJ .EQ. 0) THEN IWRK(IPIV+IWRK(IROW+I)) = IWRK(IPIV+IWRK(IROW+I)) - 1 ENDIF ENDIF 20 CONTINUE ENDIF 30 CONTINUE RETURN END * ******************************************************************** * SUBROUTINE GJWRT (K,IWRK,NWRK,IHTBL,NSTOR) * * THIS ROUTINE WRITES PIVOTING AND MATRIX REDUCTION INFORMATION * TO A FILE AS A CHECK PROCEEDURE. * * GLOBAL VARIABLES: * * 1. K - INTEGER - CURRENT STAGE IN G.E. * 2. IWRK - NWRK INTEGER ARRAY - SEE DRIVER FOR SPECIFICATION * 3. IHTBL - NSTOR INTEGER ARRAY - SEE DRIVER FOR SPECIFICATION * * GLOBAL COMMONS USED: * * 1. KUN13 - OUTPUT UNIT NUMBER * 2. IMAT - MATRIX PRINT FLAG * ******************************************************************** * INTEGER IHTBL(2,0:NSTOR),IWRK(NWRK) COMMON /WPARM/ MINFO,IHPRM,IROW,ICOL,IDEP,IPIV,LST,IPRM,IPIN, * IAPRM,KODE,IRHS,JWRK,KNUM,KGCD,IMAT,KUN13,IBASE, * ITER,ICPRM,IDIAG,ITRM,IBASIC,ITRI,INULL,IRANK,IZERO,LMAX, * NPRME,MHZ,ISOLH,KUN12 *------------------------------------------------------------------- * PRINT OUT THE PIVOT SEQUENCE *------------------------------------------------------------------- WRITE (KUN13,1050) K WRITE (KUN13,1040) IWRK(ICOL+K) WRITE (KUN13,1000) WRITE (KUN13,*) (IWRK(IPIV+IWRK(ICOL+J)),J=1,IWRK(MINFO+2)) WRITE (KUN13,1020) WRITE (KUN13,*) (IWRK(ICOL+J),J=1,IWRK(MINFO+2)) WRITE (KUN13,1030) WRITE (KUN13,*) (IWRK(IROW+J),J=1,IWRK(MINFO+1)) *------------------------------------------------------------------- * PRINT OUT THE MATRIX *------------------------------------------------------------------- IF ( IMAT .GE. 3 ) THEN WRITE (KUN13,1010) DO 20, I = 1, IWRK(MINFO+1) DO 10, J = 1, IWRK(MINFO+2)+IWRK(MINFO+3) KEY = IRHASH(IWRK(IROW+I),IWRK(ICOL+J),0,IWRK(IHPRM+1), * IHTBL,NSTOR) IF(IHTBL(2,KEY) .NE. 0 ) * WRITE (KUN13,1100) I,J,IHTBL(2,KEY) 10 CONTINUE 20 CONTINUE ENDIF * 1000 FORMAT (/' pivoting counter vector:') 1010 FORMAT (/' matrix at this stage: ') 1020 FORMAT (/' column ordering:') 1030 FORMAT (/' row ordering:') 1040 FORMAT (/' current pivot column:',I5) 1050 FORMAT (//' stage #: ',I5) 1100 FORMAT (1X,I6,I6,I11) END * ******************************************************************** * INTEGER FUNCTION INV(IX,IP) * * EUCLID'S EXTENDED ALGORITHM IS USED TO FIND THE MULTIPLICATIVE * INVERSE, INV, MODULO IP OF IX. * * *****WARNING***** * * THIS FUNCTION ASSUMES THAT IX AND IP ARE RELATIVELY PRIME. * ******************************************************************** * INTEGER IP,IR0,IX0,IR1,IX1,IR2,IX2,IQ IR0 = IP IX0 = 0 IR1 = IX IX1 = 1 10 CONTINUE IF (IABS(IR1) .EQ. 1) THEN GO TO 20 ELSE IQ = IR0/IR1 IR2 = IR0 - IQ*IR1 IX2 = IX0 - IQ*IX1 IR0 = IR1 IX0 = IX1 IR1 = IR2 IX1 = IX2 ENDIF GO TO 10 20 CONTINUE INV = MOD(IX1,IP) IF (IR1 .LT. 0) INV = -INV RETURN END * ******************************************************************** * INTEGER FUNCTION IPR (IPARM,IW1,IW2,N) * * COMPUTES AN OVER-ESTIMATE OF THE NUMBER OF PRIMES REQUIRED * FOR RSDANA. IT IS NOT NECESSARY TO CALL THIS ROUTINE IF * MPCON HAS BEEN PREVIOUSLY CALLED. * * THE ESTIMATE MAY BE GROSSLY LARGE. IT ESTIMATES THE HADAMARD * BOUND AS * * DET(A) < MAX(AIJ)*(N**(N/2)). * * INPUT GLOBAL VARIABLES: * * IPARM - (0:15) ELEMENT INTEGER ARRAY - SAME INPUT AS RSDANA. * IW1-2 - N INTEGER ARRAYS - WORK SPACE BIG ENOUGH TO STORE THE * LARGEST ELEMENT OF THE MATRIX * ******************************************************************** * INTEGER IPARM(0:15),IW1(N),IW2(N),MINFO(4) *------------------------------------------------------------------- * ESTIMATE THE PRIMES *------------------------------------------------------------------- IW2(1) = 0 REWIND (IPARM(3)) READ (IPARM(3),*) (MINFO(I),I=1,4),JBASE DO 10, K = 1, MINFO(4) READ (IPARM(3),*) I,J,NTRM READ (IPARM(3),*) (IW1(II),II=1,NTRM) IW1(1) = IW1(1)*NTRM CALL MPMAX (IW2,IW1,N,IN) IF (IN.EQ.2) CALL MPCOPY(IW2,IW1,N) 10 CONTINUE REWIND (IPARM(3)) *------------------------------------------------------------------- * ESTIMATE THE PRIMES *------------------------------------------------------------------- AMAX = 2.*FLOAT(MAX0(MINFO(1),MINFO(2))) AMAX = (AMAX/2.)*(ALOG10(AMAX)) BMAX = ALOG10(FLOAT(IW2(2))) + (ALOG10(FLOAT(IPARM(9))) * * (FLOAT(IABS(IW2(1))) - 2)) AMAX = AMAX*(1. + BMAX) BMAX = ALOG10(SQRT(FLOAT(IPARM(11))) * 0.7) IPR = INT(1.1*AMAX / BMAX) + 1 RETURN END * ******************************************************************** * INTEGER FUNCTION IRHASH (I,J,K,IHPRM,IHTBL,NSTOR) * * RECOVERS (I,J,K) DATA LOCATION IN HASHING TABLE. NOTE, K SHOULD * BE 0 WHEN THE ROW REDUCTION IS PERFORMED, AND VARIABLE ELSE. * * GLOBAL VARIABLES: * * 1. I,J,K - ARRAY LOCATION TO BE RECOVERED - WHERE * I - ROW INDEX * J - COLUMN INDEX * K - ITERATION NUMBER * 2. IHPRM - INTEGER ARRAY - CONTAINS HASHING PARAMETERS * (1) = LENGTH OF THE HASH TABLE * (2) = HASHING MULTIPLIER * (3) = ELIMINATION HASH TABLE SIZE * (4) = UNUSED HERE * (5) = COLUMNS OF AUGMENTED MATRIX * (6) = MAX. ROWS OR (COLUMNS+1) * (7) = ELIMINATION MAX. DOUBLE HASHES TRIED * (8) = UNUSED HERE * (9) = TOTAL DOUBLE HASHES TRIED * (10) = TOTAL HASH STORAGES TRIED * (11) = TERMINAL OUTPUT UNIT NUMBER * 3. IHTBL - NSTOR INTEGER ARRAY - SEE DRIVER * * LOCAL VARIABLES: * * 1. IFKY - INTEGER - KEY FUNCTION VALUE * 2. IHKY - INTEGER - HASH FUNCTION VALUE * 3. JKEY - INTEGER - KEY AT IHKY * 4. NUMH - INTEGER - ITERATION COUNTER * 5. INC - INTEGER - DOUBLE HASH INCREMENT * ******************************************************************** * INTEGER IFKY,IHKY,JKEY,INC INTEGER IHTBL(2,0:NSTOR),IHPRM(11) *------------------------------------------------------------------- * COMPUTE KEY AND HASH LOCATION - RECOVERY IS SUCESSFUL IF * HASH IS 0 OR IF HASH = KEY - ELSE DOUBLE HASH AND TRY AGAIN *------------------------------------------------------------------- IFKY = ((((K*IHPRM(6)) + (I-1))) * IHPRM(5)) + J IHKY = MOD(IFKY+(I*IHPRM(2)+J)+(I*J*K),IHPRM(3)) INC = MOD(IFKY,IHPRM(3)-2)+1 NUMH = 0 10 CONTINUE JKEY = IHTBL(1,IHKY) IF (JKEY .EQ. 0 .OR. JKEY.EQ.IFKY ) THEN GO TO 20 ELSE IHKY = MOD((IHKY+INC),IHPRM(3)) NUMH = NUMH + 1 IF (NUMH .GT. IHPRM(7)) THEN IHKY = NSTOR GO TO 20 ENDIF ENDIF GO TO 10 *------------------------------------------------------------------- * SUCCESS - ASSIGN IRHASH AND RETURN *------------------------------------------------------------------- 20 CONTINUE IRHASH = IHKY RETURN END * ******************************************************************** * INTEGER FUNCTION ISHASH (I,J,K,IHPRM,IHTBL,NSTOR) * * GETS INDEX TO STORE DATA IN THE HASHING TABLE. NOTE, K SHOULD * BE 0 WHEN THE ROW REDUCTION IS PERFORMED, AND VARIABLE ELSE. * * GLOBAL VARIABLES: * * 1. I,J,K - ARRAY LOCATION TO BE RECOVERED - WHERE * I - ROW INDEX * J - COLUMN INDEX * K - ITERATION NUMBER * 2. IHPRM - INTEGER ARRAY - CONTAINS HASHING PARAMETERS * (1) = LENGTH OF THE HASH TABLE * (2) = HASHING MULTIPLIER * (3) = ELIMINATION HASH TABLE SIZE * (4) = UNUSED HERE * (5) = COLUMNS OF AUGMENTED MATRIX * (6) = MAX. ROWS OR (COLUMNS+1) * (7) = ELIMINATION MAX. DOUBLE HASHES TRIED * (8) = UNUSED HERE * (9) = TOTAL DOUBLE HASHES TRIED * (10) = TOTAL HASH STORAGES TRIED * (11) = TERMINAL OUTPUT UNIT NUMBER * 3. IHTBL - NSTOR INTEGER ARRAY - SEE DRIVER * * LOCAL VARIABLES: * * 1. IFKY - INTEGER - KEY FUNCTION VALUE * 2. IHKY - INTEGER - HASH FUNCTION VALUE * 3. JKEY - INTEGER - KEY AT IHKY * 4. NUMH - INTEGER - ITERATION COUNTER * 5. INC - INTEGER - DOUBLE HASH INCREMENT * * ******************************************************************** * INTEGER IHTBL(2,0:NSTOR),IHPRM(11) INTEGER IFKY,IHKY,JKEY,INC *-------------------------------------------------------------------- * COMPUTE KEY AND HASH LOCATION - ALLOW STORAGE IF * 1. KEY = 0 (NEW ENTRY) * 2. THE KEY MATCHES (OVERWRITE) * 3. TABLE AT THE IHKY IS = 0 BY ROW REDUCTION * ELSE DOUBLE HASH AND TRY AGAIN *-------------------------------------------------------------------- NUMH = 0 IFKY = ((((K*IHPRM(6)) + (I-1))) * IHPRM(5)) + J IHKY = MOD(IFKY+(I*IHPRM(2)+J)+(I*J*K),IHPRM(3)) INC = MOD(IFKY,IHPRM(3)-2)+1 10 CONTINUE JKEY = IHTBL(1,IHKY) IF (JKEY.EQ.0 .OR. JKEY.EQ.IFKY .OR. IHTBL(2,IHKY).EQ.0) THEN GO TO 20 ELSE IHKY = MOD((IHKY+INC),IHPRM(3)) NUMH = NUMH + 1 IF (NUMH.GT.IHPRM(3) .AND. IHPRM(11).GT.0) THEN ITRM = IHPRM(11) WRITE (ITRM,1000) STOP ELSEIF (NUMH.GT.IHPRM(3)) THEN IHPRM(11) = -1 ENDIF GO TO 10 ENDIF *------------------------------------------------------------------- * SUCCESS - ASSIGN KEY AND RETURN - CHECK FOR MAX ITERATE *------------------------------------------------------------------- 20 CONTINUE IF (NUMH .GT. IHPRM(7)) IHPRM(7) = NUMH IF (K .EQ. 0) IHPRM(10) = IHPRM(10) + 1 IF (K .EQ. 0) IHPRM(9) = IHPRM(9) + NUMH IHTBL(1,IHKY) = IFKY ISHASH = IHKY RETURN 1000 FORMAT (//1X,'***********************************************' * ,/' MAX ITERATES EXCEEDED - HASH STRUCTURE FULL' * ,/' PARAMETER NSTOR MUST BE INCREASED TO CONTINUE' * ,/' ***********************************************') END * ******************************************************************** * INTEGER FUNCTION IWORK (IPARM) * * COMPUTES THE AMOUNT OF WORK SPACE REQUIRED FOR RSDANA. * * INPUT GLOBAL VARIABLES: * * IPARM - (0:15) ELEMENT INTEGER ARRAY - SAME INPUT AS RSDANA. * ******************************************************************** * INTEGER IPARM(0:15),LMAX,MAXP * MAXP = INT(SQRT(FLOAT(IPARM(11)))) LMAX = 2*(INT( FLOAT(IPARM(12)) * ALOG10(FLOAT(MAXP)/(2.0)) / * ALOG10(FLOAT(IPARM(9))) ) + 10) * IWORK = 1000 + (3*IPARM(0)) + (2*IPARM(1)) + (3*IPARM(12)) + * (2*MAX0(IPARM(0),IPARM(1)+IPARM(2))) + (12*LMAX) * RETURN END * ******************************************************************** * SUBROUTINE JRHASH (I,J,L,ITN,IOV,JZ,IWRK,NWRK,IHTBL,NSTOR) * * RECOVERS (I,J) SOLUTION DATA ELEMENT IN SOLUTION HASHING TABLE * RETURNS IT AT L POSITION IN IWRK. REQUIRED INPUT IS A * NON-ZERO VALUE I AND J, AND SOME VALUE OF ITN. * * GLOBAL VARIABLES: * * 1. I,J - ROW AND COLUMN SOLUTION ELEMENT TO RECOVER * 2. L - ADDRESS OF THE OUTPUT IN IWRK * 3. ITN - MAXIMUM NUMBER OF TERMS * 4. IOV - FLAG TO ZERO THE ENTRY ONCE RECOVERED * = 0 THEN ZERO * = ANYTHING ELSE THEN DONT * 5. JZ - SIGNALS WHEN ELEMENT RECOVERED IS ZERO * = 0 IF (L = 0) * = 1 ELSE * * LOCAL VARIABLES: * * 1. IFKY - KEY FUNCTION VALUE * 2. IHKY - HASH FUNCTION VALUE * 3. JKEY - KEY FUNCTION VALUE * 4. INC - DOUBLE HASH INCREMENT * 5. NUMH - DOUBLE HASH COUNTER * * ******************************************************************** * INTEGER IHTBL(2,0:NSTOR),IWRK(NWRK),L,ITN,IOV,JZ INTEGER IFKY,IHKY,JKEY,INC,NUMH,LK,LK1 COMMON /WPARM/ MINFO,IHPRM,IROW,ICOL,IDEP,IPIV,LST,IPRM,IPIN, * IAPRM,KODE,IRHS,JWRK,KNUM,KGCD,IMAT,KUN13,IBASE, * ITER,ICPRM,IDIAG,ITRM,IBASIC,ITRI,INULL,IRANK,IZERO,LMAX, * NPRME,MHZ,ISOLH,KUN12 *------------------------------------------------------------------- * INITIALIZE - COMPUTE IFKY, MULTIPLIER, IHKY AND DOUBLE HASH * INCREMENT *------------------------------------------------------------------- JZ = 0 IF ( ITN .EQ. 0 ) RETURN JKEY = (I-1)*IWRK(IHPRM+5) + J DO 20, K = 1, (ITN+MOD(ITN,2))/2 NUMH = 0 LK = 2*K LK1 = LK - 1 IFKY = JKEY + (K*IWRK(IHPRM+5)*IWRK(MINFO+1)) IHKY = MOD(K*((I*IWRK(IHPRM+2))+J),IWRK(IHPRM+4)) + MHZ INC = MOD(IFKY,IWRK(IHPRM+4)-2)+1 *------------------------------------------------------------------- * FIND ZERO TERMS *------------------------------------------------------------------- 10 CONTINUE IF ( IHTBL(1,IHKY).EQ.0 ) THEN IWRK(L+LK1) = 0 IWRK(L+LK) = 0 *------------------------------------------------------------------- * SUCCESSFUL HASH WITH NONZERO TERMS - UNPACK THE TERMS *------------------------------------------------------------------- ELSEIF ( IHTBL(1,IHKY).EQ.IFKY ) THEN JZ = 1 KK = IHTBL(2,IHKY) IWRK(L+LK) = MOD(KK,IWRK(IPRM+1)) IWRK(L+LK1) = (KK-IWRK(L+LK)) / IWRK(IPRM+1) IF ( IWRK(L+LK) .GE. (IWRK(IAPRM+LK)+1)/2 ) * IWRK(L+LK) = IWRK(L+LK)-IWRK(IAPRM+LK) IF ( IWRK(L+LK1) .GE. (IWRK(IAPRM+LK1)+1)/2 ) * IWRK(L+LK1) = IWRK(L+LK1)-IWRK(IAPRM+LK1) IF ( IOV .EQ. 0 ) THEN IHTBL(1,IHKY) = -1 IHTBL(2,IHKY) = 0 ENDIF *------------------------------------------------------------------- * UNSUCCESSFUL HASH - TRY AGAIN IF LIMIT NOT EXCEEDED *------------------------------------------------------------------- ELSE IHKY = MOD((IHKY-MHZ)+INC,IWRK(IHPRM+4)) + MHZ NUMH = NUMH + 1 IF ( NUMH .GT. IWRK(IHPRM+8) ) THEN IWRK(L+LK) = 0 IWRK(L+LK1) = 0 GO TO 20 ENDIF GO TO 10 ENDIF 20 CONTINUE RETURN END * ******************************************************************** * SUBROUTINE JSHASH (I,J,L,ITN,IWRK,NWRK,IHTBL,NSTOR) * * STORES THE SOLUTION ELEMENT IN IWRK(L) - IWRK(L+ITN) * IN THE SOLUTION HASHING TABLE. REQUIRED INPUT IS A NONZERO * VALUE I AND J, AND SOME NON-NEGATIVE VALUE FOR ITN. * * GLOBAL VARIABLES: * * 1. I,J - ROW AND COLUMN MATRIX LOCATION TO STORE * 2. L - ADDRESS OF ELEMENT IN IWRK * 3. ITN - NUMBER OF TERMS OF L TO STORE * 4. IWRK - NWRK INTEGER ARRAY - SEE DRIVER * 5. IHTBL - NSTOR INTEGER ARRAY - SEE DRIVER * * LOCAL VARIABLES: * * 1. IFKY - KEY FUNCTION VALUE * 2. IHKY - HASH FUNCTION VALUE * 3. JKEY - KEY FUNCTION VALUE * 4. INC - DOUBLE HASH INCREMENT * 5. NUMH - ITERATION COUNTER * * ******************************************************************** * INTEGER IHTBL(2,0:NSTOR),IWRK(NWRK),I,J,L,ITN INTEGER IFKY,IHKY,JKEY,INC,LK,LK1 COMMON /WPARM/ MINFO,IHPRM,IROW,ICOL,IDEP,IPIV,LST,IPRM,IPIN, * IAPRM,KODE,IRHS,JWRK,KNUM,KGCD,IMAT,KUN13,IBASE, * ITER,ICPRM,IDIAG,ITRM,IBASIC,ITRI,INULL,IRANK,IZERO,LMAX, * NPRME,MHZ,ISOLH,KUN12 *------------------------------------------------------------------- * INITIALIZE - COMPUTE IFKY, MULTIPLIER, IHKY AND DOUBLE HASH * INCREMENT - SKIP THE ZERO TERMS *------------------------------------------------------------------- JKEY = (I-1)*IWRK(IHPRM+5) + J IWRK(L+ITN+1) = 0 DO 20, K = 1, (ITN+MOD(ITN,2))/2 NUMH = 0 LK = 2*K LK1 = LK - 1 IF ( (IWRK(L+LK).EQ.0) .AND. (IWRK(L+LK1).EQ.0) ) GO TO 20 IFKY = JKEY + (K*IWRK(IHPRM+5)*IWRK(MINFO+1)) IHKY = MOD(K*((I*IWRK(IHPRM+2))+J),IWRK(IHPRM+4)) + MHZ INC = MOD(IFKY,IWRK(IHPRM+4)-2)+1 *-------------------------------------------------------------------- * STORE TERMS - ALLOW STORAGE IF KEY <= 0 OR IF THE KEY MATCHES. * WHEN STORAGE IS ALLOWED PACK THE TERMS SO THAT TWO CONSECUTIVE * TERM OCCUPY ONE HASH LOCATION. *-------------------------------------------------------------------- 10 CONTINUE IF (IHTBL(1,IHKY).LE.0 .OR. IHTBL(1,IHKY).EQ.IFKY) THEN IHTBL(1,IHKY) = IFKY ILK = IWRK(L+LK) ILK1 = IWRK(L+LK1) IF ( ILK .LT. 0 ) ILK = ILK + IWRK(IAPRM+LK) IF ( ILK1 .LT. 0 ) ILK1 = ILK1 + IWRK(IAPRM+LK1) IHTBL(2,IHKY) = (ILK1*IWRK(IPRM+1)) + ILK *------------------------------------------------------------------- * DOUBLE HASH UNTIL SUCCESS - CHECK TO INSURE THE SOLUTION HASH * TABLE IS NOT FULL *------------------------------------------------------------------- ELSE IHKY = MOD((IHKY-MHZ)+INC,IWRK(IHPRM+4)) + MHZ NUMH = NUMH + 1 IF (NUMH .GT. IWRK(IHPRM+4)) THEN ISOLH = 1 IWRK(IHPRM+8) = IWRK(IHPRM+4) RETURN ELSEIF (NUMH .GT. IWRK(IHPRM+8)) THEN IWRK(IHPRM+8) = NUMH ENDIF GO TO 10 ENDIF 20 CONTINUE RETURN END * SUBROUTINE MPADD (IX,IY,IZ,N,BASE) * * TO ADD MP NON NEGATIVE INTEGERS IX + IY = IZ. * * GLOBAL VARIABLES: * * 1. IX - N INTEGER ARRAY - COEFFIENTS OF IX * 2. IY - N INTEGER ARRAY - COEFFIENTS OF IY * 3. IZ - N INTEGER ARRAY - COEFFIENTS OF IZ * 4. N - INTEGER - NUMBER OF ARRAY ELEMENTS * 5. BASE - INTEGER - NUMBER BASE OF IX,IY,IZ * * REF: D.E. KNUTH, 'THE ART OF COMPUTER PROGRAMMING', * SEC. 4.3.1 * * ******************************************************************** * INTEGER IX(N),IY(N),IZ(N),BASE * NX = IABS(IX(1)) NY = IABS(IY(1)) NZ = MAX0(NX,NY) MZ = MIN0(NX,NY) K = 0 DO 10, J = 0, MZ-2 INT = IX(NX-J) + IY(NY-J) + K IF (INT .LT. BASE) THEN IZ(NZ-J) = INT K = 0 ELSE IZ(NZ-J) = INT - BASE K = 1 ENDIF 10 CONTINUE IF ( NZ .EQ. MZ ) THEN IZ(2) = (K*BASE) + IZ(2) ELSEIF ( NZ .EQ. NX ) THEN IZ(NZ-MZ+1) = IX(NZ-MZ+1) + K DO 20, J= NZ-MZ, 2, -1 IZ(J) = IX(J) 20 CONTINUE ELSE IZ(NZ-MZ+1) = IY(NZ-MZ+1) + K DO 30, J= NZ-MZ, 2, -1 IZ(J) = IY(J) 30 CONTINUE ENDIF IZ(1) = NZ IF (IZ(2) .EQ. 0) CALL MPPACK(IZ,N) RETURN END * ******************************************************************** * SUBROUTINE MPBASE(IX,N,MBASE,NBASE,IW1) * * TO CONVERT MULTIPLE PRECISION IX RADIX MBASE TO * IX RADIX NBASE. MBASE AND NBASE MUST BE ONE WORD * INTEGERS, SO THAT THEIR SQUARE IS LESS THAN THE * COMPUTER WORD. * * GLOBAL VARIABLES: * * 1. IX - N INTEGER ARRAY - NUMBER TO BE CONVERTED ON * INPUT - SAME NUMBER DIFFERENT BASE ON OUTPUT. * 2. N - INTEGER - NUMBER OF ARRAY ELEMENTS IN IX * 3. MBASE - BASE ON INPUT * 4. NBASE - BASE ON OUTPUT. * ******************************************************************** * INTEGER IX(N),MBASE,NBASE,IW1(N) ISIGN = 1 IF ( IX(1) .LT. 0 ) ISIGN = -1 *-------------------------------------------------------------------- * PERFORM DIVISION OF IX WITH NBASE - RESULTS IN QUOTIENT AND * REMAINDER = JTH DECIMAL OF IX BASE(NBASE) *-------------------------------------------------------------------- J = 0 NTRM = IABS(IX(1)) 10 CONTINUE INT = IX(2) J = J+1 DO 20, I = 2,NTRM,1 IREM = MOD(INT,NBASE) IX(I) = (IX(I) - IREM) / NBASE IX(I+1) = IX(I+1) + (IREM*MBASE) INT = IX(I+1) 20 CONTINUE IW1(J) = IREM *-------------------------------------------------------------------- * DONE IF QUOTIENT IS ZERO *-------------------------------------------------------------------- IMAX = 0 DO 30, I = 2, NTRM IMAX = MAX0(IMAX,IX(I)) 30 CONTINUE IF (IMAX .GT. 0) GO TO 10 *-------------------------------------------------------------------- * EXCHANGE TEMP WITH IX - RETURN *-------------------------------------------------------------------- DO 40 I = 1, J IX(J-I+2) = IW1(I) 40 CONTINUE IX(1) = ISIGN * (J+1) IF ( IX(2) .EQ. 0 ) CALL MPPACK (IX,N) RETURN END * ******************************************************************** * SUBROUTINE MPCON (MPARM,IWRK,NWRK) * * READS A REGULAR FORMAT OR ARBITRARY PRECISION INPUT DATA * FILE WITH RATIONAL COEFFICIENTS OF A MATRIX. CONVERTS * THE RATIONALS TO INTEGERS ROW-WISE. IT THEN OUTPUTS TO A * FILE THAT IS USED AS INPUT TO RSDANA. * * ****INPUTS * * GLOBAL VARIABLES: * * 1. MPARM - 7 INTEGER ARRAY - CONTROL PARAMETERS WHERE * (1) - INPUT FILE UNIT NUMBER * (2) - OUTPUT FILE UNIT NUMBER * (3) - INTEGER - OUTPUT RADIX * (4) - FLAG TO SIGNAL HOW TO READ INPUT * = 0 => ARBITRARY PRECISION * = 1 => REGULAR FORMAT * (5) - TERMINAL OUTPUT FLAG * = 0 => NO INTERMEDIATE OUTPUT * > 1 => OUTPUT * (6) - TERMINAL OUTPUT UNIT NUMBER * (7) - COMPUTER WORD * * 2. IWRK - NWRK INTEGER ARRAY - LARGEST WORK SPACE * AVAILABLE (USUALLY IT WILL BE THE HASH TABLE) * * INPUTS FILE STRUCTURE: EITHER DATA FILE TYPE SHOULD CONTAIN * THE AUGMENTED MATRIX SORTED BY ROWS INCLUDING THE RIGHT HAND * SIDES NUMBERED AS THEY WOULD BE IN THE AUGMENTED MATRIX. * ZERO ELEMENTS SHOULD NOT BE INCLUDED. THE ALLOWED DATA * STURCTURES ARE; * * * 1. REGULAR DATA FILE: * * NROW, NCOL, NRHS, NELE * I, J, NIJ, DIJ - REPEATED NELE TIMES * * * 2. ARBITRARY PRECISION DATA FILE: * * NROW, NCOL, NRHS, NELE * I, J, ] * NIJ ] -- REPEATED NELE TIMES * DIJ ] * * WHERE: NROW = # OF ROWS * NCOL = # OF COLUMNS * NRHS = # OF RIGHT HAND SIDES * NELE = # OF NON-ZERO ELEMENTS * (I,J) = IS THE MATRIX INDEX * NIJ = THE NUMERATOR OF THE (I,J) ELEMENT * DIJ = THE DEMONINATOR OF THE (I,J) ELEMENT * * ***************************************************************** * * VERY IMPORTANT FOR ARBIRARY PRECISON DATA: EVERY NIJ, DIJ * * * MUST BE TERMINATED WITH A "$" * * ***************************************************************** * * * EXAMPLE OF A REGULAR DATA FILE TO SOLVE AX = B WHERE A IS A * 3X2 MATRIX AND B IS A 3X1 VECTOR * * 5*A = [ 1 0 ] * [ 0 3 ] AND 5*B = [ 1, 0, 5 ]' * [ 0 5 ] * * THEN THE DATA FILE SHOULD CONTAIN * * 3, 2, 1, 5 (THE MATRIX INFORMATION) * 1, 1, 1, 5 (I.E. THE 1,1 ELEMENT OF A) * 1, 3, 1, 5 (I.E. THE 1,1 ELEMENT OF B OR 1,3 ELEMENT OF A AUG.) * 2, 2, 3, 5 ( ... ) * 3, 2, 1, 1 * 3, 3, 1, 1 * * ****OUTPUT * * * 1. IF AN ERROR HAS OCCURED IN INPUT MPARM(1) WILL BE SET * TO -1 * * 2. MPARM(7) IS OVERWRITTEN WITH AN OVER ESTIMATE * OF THE NUMBER OF PRIMES RSDANA WILL USE TO FIND * THE EXACT SOLUTION OF THE PROBLEM. UNFORTUNATELY * THE ESTIMATE MAY BE GROSSLY LARGE. IT ESTIMATES * THE HADAMARD BOUND AS * * DET(A) < MAX(AIJ)*(N**(N/2)). * * 3. OUTPUT DATA FILE - STRUCTURE CONTAINS * * NROW, NCOL, NRHS, NELE, JBASE * I, J, NTRM ] * (AIJ)K K=1,2,...,NTRM ] - REPEATED NELE TIMES * 0, 0, NTRM ] * (D)K K=1,2,...,NTRM ] - DIAGONAL DETERMINANT * * WHERE: NROW = # OF ROWS * NCOL = # OF COLUMNS * NRHS = # OF RIGHT HAND SIDES * NELE = # OF NON-ZERO ELEMENTS * I,J = ELEMENT MATRIX INDEX * NTRM = NUMBER OF TERMS RADIX(JBASE) OF AIJ * (AIJ)K = AIJ ELEMENT WHERE K IS THE TERM INDEX * AND (AIJ)1 REPRESENTS THE SIGN OF THE ELEMENT. * (D)K = DETERMINANT OF THE DIAGONAL MATRIX WHICH * TRANSFORMS A INTO AN INTEGER MATRIX. STRUCTURE * OF (D)K IS THE SAME AS (AIJ)K. * * THE RECONSTRUCTION OF AN MP NUMBER IS * * AIJ = (AIJ)1 * ( (AIJ)2*(JBASE**(NTMR-2)) + * (AIJ)3*(JBASE**(NTMR-3)) + * ... * (AIJ)NTRM-1*JBASE + (AIJ)NTRM ) * * * TO READ THE DATA INTO A 3-DIMENSIONAL ARRAY AND CONVERT IT * TO MP FORMAT FOLLOW THE SAMPLE OF CODE * * READ (I,*) NROW,NCOL,NRHS,NELE,JBASE * DO 10, L = 1, NELE * READ (I,*) I,J,NTRM * READ (I,*) (A(I,J,K),K=1,NTRM) * A(I,J,1) = A(I,J,1)*NTRM * 10 CONTINUE * ******************************************************************** * INTEGER IWRK(NWRK),MPARM(7) INTEGER MINFO(4) LOGICAL TEST CHARACTER*1 A(80) *------------------------------------------------------------------- * INPUT PARAMETER CHECKS - I/O UNIT NUMBERS *------------------------------------------------------------------- IER = 0 TEST = .TRUE. IF (MPARM(6) .LE. 0) THEN MPARM(6) = 3 WRITE (MPARM(6),1100) WRITE (MPARM(6),1110) IER = 1 ENDIF IF (MPARM(1).LE.0) THEN TEST = .FALSE. ELSE INQUIRE (MPARM(1),OPENED=TEST) ENDIF IF ( TEST ) GO TO 10 IF ( IER .EQ. 0 ) WRITE (MPARM(6),1100) WRITE (MPARM(6),1120) IER = 1 10 CONTINUE IF (MPARM(2).LE.0) THEN TEST = .FALSE. ELSE INQUIRE (MPARM(2),OPENED=TEST) ENDIF IF ( TEST ) GO TO 20 IF ( IER .EQ. 0 ) WRITE (MPARM(6),1100) WRITE (MPARM(6),1130) IER = 1 20 CONTINUE *------------------------------------------------------------------- * INPUT PARAMETER CHECKS - I/O UNIT NUMBERS *------------------------------------------------------------------- IF ( MPARM(3).LT.2 ) THEN IF ( IER .EQ. 0 ) WRITE (MPARM(6),1100) WRITE (MPARM(6),1140) IER = 1 ENDIF IF (MPARM(4).LT.0 .OR. MPARM(4).GT.1) THEN IF ( IER .EQ. 0 ) WRITE (MPARM(6),1100) WRITE (MPARM(6),1150) IER = 1 ENDIF IF (MPARM(5).LT.0 ) THEN IF ( IER .EQ. 0 ) WRITE (MPARM(6),1100) WRITE (MPARM(6),1160) IER = 1 ENDIF IF ( IER .EQ. 1 ) THEN MPARM(1) = -1 WRITE (MPARM(6),1170) RETURN ENDIF *------------------------------------------------------------------- * INITIALIZE WORK SPACE *------------------------------------------------------------------- DO 30, I = 1,NWRK IWRK(I) = 0 30 CONTINUE *------------------------------------------------------------------- * READ THE MATRIX PARAMETERS - FIND THE MAX NUMBER OF NONZERO * ELEMENTS IN ANY COLUMN. *------------------------------------------------------------------- READ (MPARM(1),*) (MINFO(I), I=1,4) WRITE (MPARM(2),*) (MINFO(I), I=1,4),MPARM(3) IF (MINFO(1) .EQ. 1) THEN MCOL = MINFO(2) + 2 ELSE MCOL = 3 ENDIF IR = 0 JC = 0 DO 70, I = 1, MINFO(4) IF ( MPARM(4) .EQ. 1 ) THEN READ (MPARM(1),*) II,J,KK,LL ELSE READ (MPARM(1),*) II,J *------------------------------------------------------------------- * NEED TO FIND THE NEXT PAIR I,J *-------------------------------------------------------------------- DO 60, KTEMP =1,2 40 CONTINUE READ (MPARM(1),1000) (A(JTEMP),JTEMP=1,80) DO 50, JTEMP = 1,80 IF (A(JTEMP).EQ.'$') GO TO 60 50 CONTINUE GO TO 40 60 CONTINUE ENDIF IF (II .EQ. IR) THEN JC = JC + 1 ELSE MCOL = MAX(MCOL,JC) JC = 1 IR = II ENDIF 70 CONTINUE MCOL = MCOL + 1 *------------------------------------------------------------------- * SPLIT THE WORK SPACE *------------------------------------------------------------------- NVEC = (2*MCOL) + 15 NTRM = NWRK / NVEC LDNM = (MCOL * NTRM) ICOL = (2 * MCOL * NTRM) + 1 LCMD = ICOL + NTRM LST = LCMD + NTRM IW1 = LST + NTRM IW2 = IW1 + NTRM IW3 = IW2 + NTRM IW4 = IW3 + NTRM IW5 = IW4 + NTRM IW6 = IW5 + NTRM IDET = IW6 + NTRM + NTRM IEST = IDET + NTRM + NTRM IWRK(IDET) = 2 IWRK(IDET+1) = 1 *------------------------------------------------------------------- * READ MATRIX ONE ROW AT A TIME - COMPUTE LCM OF THE * DENOMINATORS PAIRWISE. *------------------------------------------------------------------- REWIND (MPARM(1)) READ (MPARM(1),*) (MINFO(I), I=1,4) L = 0 IROW = 1 IF (MPARM(5).GE.1) WRITE (MPARM(6),1020) * IF (MPARM(5).GE.1) WRITE (MPARM(6),1030) IROW IWRK(LCMD) = 2 IWRK(LCMD+1) = 1 DO 120 I = 1, MINFO(4) IF ( MPARM(4) .EQ. 0 ) THEN READ (MPARM(1),*) IR,IC ELSE READ (MPARM(1),*) IR,IC,INUM,IDEN ENDIF IF ( IR.EQ.IROW ) THEN L = L+1 IWRK(ICOL+L) = IC NUMER = ((L-1)*NTRM) + 1 IDENM = ((L-1)*NTRM) + 1 + LDNM IF ( MPARM(4) .EQ. 0 ) THEN CALL MPREAD (IWRK(NUMER),NTRM,MPARM(1),MPARM(3), * IWRK(LST)) CALL MPREAD (IWRK(IDENM),NTRM,MPARM(1),MPARM(3), * IWRK(LST)) ELSE CALL MPMAKE (INUM,IWRK(NUMER),NTRM,MPARM(3)) CALL MPMAKE (IDEN,IWRK(IDENM),NTRM,MPARM(3)) ENDIF DO 80, J = LST, IDET-1, NTRM IWRK(J) = 0 80 CONTINUE CALL MPCOPY (IWRK(LST),IWRK(IDENM),NTRM) CALL MPLCM (IWRK(LST),IWRK(LCMD),NTRM,MPARM(3),IWRK(IW1), * IWRK(IW3),IWRK(IW4),IWRK(IW5),IWRK(IW6)) *------------------------------------------------------------------- * NEW ROW IS FOUND - MULTIPLY OLD ROW BY LCMD AND OUTPUT *------------------------------------------------------------------- ELSE DO 100 J = 1, L, 1 DO 90, JJ = LST, IDET-1, NTRM IWRK(JJ) = 0 90 CONTINUE ISIGN = 1 NUMER = ((J-1)*NTRM) + 1 IDENM = ((J-1)*NTRM) + 1 + LDNM IF ( IWRK(NUMER)*IWRK(IDENM) .LT. 0 ) ISIGN = -1 CALL MPMULT (IWRK(NUMER),IWRK(LCMD),IWRK(LST),NTRM, * MPARM(3)) CALL MPDIV (IWRK(LST),IWRK(IDENM),IWRK(NUMER),IWRK(IW2), * NTRM,MPARM(3),IWRK(IW4),IWRK(IW6)) NN = IABS(IWRK(NUMER)) WRITE (MPARM(2),*) IROW,IWRK(ICOL+J),NN WRITE (MPARM(2),1010) ISIGN,(IWRK(NUMER+K),K=1,NN-1) WRITE (MPARM(2),1020) CALL MPMAX (IWRK(NUMER),IWRK(IEST),NTRM,NN) IF (NN.EQ.1) CALL MPCOPY (IWRK(IEST),IWRK(NUMER),NTRM) 100 CONTINUE *------------------------------------------------------------------- * INITIALIZE NEW ROW *------------------------------------------------------------------- CALL MPMULT (IWRK(IDET),IWRK(LCMD),IWRK(LST),NTRM,MPARM(3)) CALL MPCOPY (IWRK(IDET),IWRK(LST),NTRM) DO 110, JJ = 1, IDET-1, NTRM IWRK(JJ) = 0 110 CONTINUE L = 1 IROW = IR IWRK(ICOL+1) = IC NUMER = 1 IDENM = LDNM + 1 * IF (MPARM(5).GE.1) WRITE (MPARM(6),1030) IROW IF ( MPARM(4) .EQ. 0 ) THEN CALL MPREAD (IWRK(NUMER),NTRM,MPARM(1),MPARM(3), * IWRK(LST)) CALL MPREAD (IWRK(IDENM),NTRM,MPARM(1),MPARM(3), * IWRK(LST)) ELSE CALL MPMAKE (INUM,IWRK(NUMER),NTRM,MPARM(3)) CALL MPMAKE (IDEN,IWRK(IDENM),NTRM,MPARM(3)) ENDIF CALL MPCOPY (IWRK(LCMD),IWRK(IDENM),NTRM) ENDIF 120 CONTINUE *------------------------------------------------------------------- * END OF FILE - MULTIPLY LAST ROW BY LCMD AND OUTPUT *------------------------------------------------------------------- DO 140 J = 1, L DO 130, JJ = LST, IDET-1, NTRM IWRK(JJ) = 0 130 CONTINUE ISIGN = 1 NUMER = ((J-1)*NTRM) + 1 IDENM = ((J-1)*NTRM) + 1 + LDNM IF ( IWRK(NUMER)*IWRK(IDENM) .LT. 0 ) ISIGN = -1 CALL MPMULT (IWRK(NUMER),IWRK(LCMD),IWRK(LST),NTRM, * MPARM(3)) CALL MPDIV (IWRK(LST),IWRK(IDENM),IWRK(NUMER),IWRK(IW2), * NTRM,MPARM(3),IWRK(IW4),IWRK(IW6)) NN = IABS(IWRK(NUMER)) WRITE (MPARM(2),*) IROW,IWRK(ICOL+J),NN WRITE (MPARM(2),1010) ISIGN,(IWRK(NUMER+K),K=1,NN-1) WRITE (MPARM(2),1020) CALL MPMAX (IWRK(NUMER),IWRK(IEST),NTRM,NN) IF (NN.EQ.1) CALL MPCOPY (IWRK(IEST),IWRK(NUMER),NTRM) 140 CONTINUE *------------------------------------------------------------------- * OUTPUT THE DENOMINATOR OF THE DIAGONAL MATRIX *------------------------------------------------------------------- CALL MPMULT (IWRK(IDET),IWRK(LCMD),IWRK(LST),NTRM,MPARM(3)) NN = IABS(IWRK(LST)) WRITE (MPARM(2),*) 0, 0, NN WRITE (MPARM(2),1010) 1,(IWRK(LST+K),K=1,NN-1) WRITE (MPARM(2),1020) WRITE (MPARM(2),1040) REWIND (MPARM(2)) *------------------------------------------------------------------- * ESTIMATE THE PRIMES *------------------------------------------------------------------- AMAX = FLOAT(MAX0(MINFO(1),MINFO(2))) AMAX = (AMAX/2.)*(ALOG10(AMAX)) BMAX = ALOG10(FLOAT(IWRK(IEST+1))) + (ALOG10(FLOAT(MPARM(3))) * * (FLOAT(IABS(IWRK(IEST))) - 2)) AMAX = AMAX*(1. + BMAX) BMAX = ALOG10(SQRT(FLOAT(MPARM(7))) * 0.7) MPARM(7) = INT(1.1*AMAX/BMAX) + 1 RETURN 1000 FORMAT (80A1) 1010 FORMAT (11(10I8/)) 1020 FORMAT (' ') 1030 FORMAT (' preprocessing row ',I5) 1040 FORMAT (' final element is determinant of scaling matrix') 1100 FORMAT (/5X, * '********************************************************',/5X, * ' ',/5X, * ' input error(s) encountered in mpcon' * ' ',/5X) 1110 FORMAT (/5X,'terminal unit number error ') 1120 FORMAT (/5X,'matrix file unit error - unit not opened') 1130 FORMAT (/5X,'output file unit error - unit not opened') 1140 FORMAT (/5X,'base error ') 1150 FORMAT (/5X,'file type must equal 0 or 1') 1160 FORMAT (/5X,'intermediate output parameter error - allowed' * ' range (0-4)') 1170 FORMAT (/5X, * '********************************************************',/) END * ******************************************************************** * SUBROUTINE MPCOPY (IX,IY,N) * * COPY CONTENTS OF IY INTO IX. * * GLOBAL VARIABLES: * 1. IX - N INTEGER ARRAY - COEFFICIENTS OF IX * 2. IY - N INTEGER ARRAY - COEFFICIENTS OF IY * 3. N - INTEGER - NUMBER OF ARRAY ELEMENTS * ******************************************************************** * INTEGER IX(N),IY(N) * IF (IABS(IY(1)) .LE. 1) THEN IX(1) = 0 ELSE DO 20, I = 1, IABS(IY(1)) IX(I) = IY(I) 20 CONTINUE ENDIF RETURN END * ******************************************************************** * SUBROUTINE MPDIV (IX,IY,IQ,IR,N,IBASE,JU,JV) * * TO DIVIDE MP INTEGERS IX / IY. RESULTS IN MP QUOTIENT, IQ, * AND MP REMAINDER, IR. * * GLOBAL VARIABLES: * * 1. IX - N INTEGER ARRAY - COEFFIENTS OF DIVIDEND * 2. IY - N INTEGER ARRAY - COEFFIENTS OF DIVISOR * 3. IQ - N INTEGER ARRAY - COEFFIENTS OF QUOTIENT * 4. IR - N INTEGER ARRAY - COEFFIENTS OF REMAINDER * 5. N - INTEGER - NUMBER OF ARRAY ELEMENTS * 6. IBASE - INTEGER - NUMBER IBASE OF IX,IY,IQ,IR * 7. JU,JV - WORK SPACE * * REF: D.E. KNUTH, 'the art of computer programming', * SEC. 4.3.1 * ******************************************************************** * INTEGER IX(N),IY(N),IQ(N),IR(N),IBASE,JU(0:N),JV(0:N) INTEGER VHAT,QHAT *------------------------------------------------------------------- * EASY CASES FIRST *------------------------------------------------------------------- CALL MPMAX (IX,IY,N,IN) IF ( IN .EQ. 0 ) THEN IF ( IX(1)*IY(1) .GT. 0 ) THEN IQ(1) = 2 ELSE IQ(1) = -2 ENDIF IQ(2) = 1 IR(1) = 0 RETURN ELSEIF ( IN .EQ. 2 ) THEN IQ(1) = 0 CALL MPCOPY (IR,IX,N) RETURN ELSEIF (IY(1).EQ.2 .AND. IY(2).EQ.1) THEN CALL MPCOPY (IQ,IX,N) IR(1) = 0 RETURN ENDIF *------------------------------------------------------------------- * INITIALIZE *------------------------------------------------------------------- ISIGN = 1 IF (IX(1)*IY(1) .LT. 0) ISIGN = -1 NX = IABS(IX(1)) NY = IABS(IY(1)) NQ = NX - NY + 1 NR = NY - 1 IQ(1) = (NQ+1)*ISIGN IR(1) = (NR+1)*ISIGN DO 10, I = 2, NQ+1 IQ(I) = 0 10 CONTINUE DO 20, I = 2, NR+1 IR(I) = 0 20 CONTINUE DO 30, I = 0, NX+2 JU(I) = 0 JV(I) = 0 30 CONTINUE *------------------------------------------------------------------- * NORMALIZE *------------------------------------------------------------------- NORM = (IBASE - MOD(IBASE,IY(2)+1)) / (IY(2) + 1) K = 0 DO 40, I = NX, 2, -1 INT = IX(I)*NORM + K JU(I-1) = MOD(INT,IBASE) K = (INT - JU(I-1)) / IBASE 40 CONTINUE JU(0) = K K = 0 DO 50, I = NY, 2, -1 INT = IY(I)*NORM + K JV(I-1) = MOD(INT,IBASE) K = (INT - JV(I-1)) / IBASE 50 CONTINUE *------------------------------------------------------------------- * DIVIDE JU(J:J+N) BY JV(1:N) FOR J = 0, NR. CALCULATE AN * ESTIMATE TO THE QUOTIENT FIRST (QHAT). *------------------------------------------------------------------- DO 100, J = 0, NQ-1 INT = (JU(J)*IBASE) + JU(J+1) IF (JU(J) .EQ. JV(1)) THEN QHAT = IBASE - 1 ELSE QHAT = (INT - MOD(INT,JV(1))) / JV(1) ENDIF 60 CONTINUE IF (JV(2)*QHAT .GT. ((INT-(QHAT*JV(1)))*IBASE) * +JU(J+2)) THEN QHAT = QHAT - 1 ELSE GO TO 70 ENDIF GO TO 60 *------------------------------------------------------------------- * MULTIPLY AND SUBTRACT *------------------------------------------------------------------- 70 CONTINUE KV = 0 KU = 0 DO 80, I = NR, 1, -1 INT = JV(I)*QHAT + KV VHAT = MOD(INT,IBASE) KV = (INT - VHAT) / IBASE INT = JU(I+J) - VHAT + KU IF (INT .LT. 0) THEN JU(I+J) = INT + IBASE KU = -1 ELSE JU(I+J) = INT KU = 0 ENDIF 80 CONTINUE JU(J) = JU(J) + KU *------------------------------------------------------------------- * TEST REMAINDER *------------------------------------------------------------------- IQ(J+2) = QHAT IF (JU(J) .LT. 0) THEN IQ(J+2) = IQ(J+2) - 1 K = 0 DO 90, I = NR, 1, -1 INT = JU(J+I) + JV(I) + K IF ( INT .LT. IBASE ) THEN JU(J+I) = INT K = 0 ELSE JU(J+I) = INT - IBASE K = 1 ENDIF 90 CONTINUE JU(J) = JU(J) + K - IBASE ENDIF 100 CONTINUE *------------------------------------------------------------------- * UNNORMALIZE AND COMPUTE THE REMAINDER *------------------------------------------------------------------- IF (IQ(2) .EQ. 0) CALL MPPACK(IQ,N) K = 0 DO 110, I = NQ,NX-1 JU(I) = (K*IBASE) + JU(I) K = MOD(JU(I),NORM) IR(I-NQ+2) = (JU(I) - K) / NORM 110 CONTINUE IF (IABS(IR(1)).GT.1 .AND. IR(2).EQ.0) CALL MPPACK(IR,N) CALL MPMAX (IR,IY,N,IN) IF ( IN .EQ. 0 ) THEN IR(1) = 0 IQ(IABS(IQ(1))) = IQ(IABS(IQ(1))) + 1 ENDIF DO 130, I = IABS(IQ(1)), 3, -1 120 CONTINUE IF ( IQ(I) .GE. IBASE ) THEN IQ(I) = IQ(I) - IBASE IQ(I-1) = IQ(I-1) + 1 G O TO 120 ENDIF 130 CONTINUE RETURN END * ******************************************************************** * SUBROUTINE MPGCD (IX,IGCD,N,IBASE,IW1,IW2,JU,JV) * * EUCLIDEAN ALGORITHM TO FIND THE GREATEST COMMON DIVISOR * OF MP INTEGERS IX, AND IGCD. THE ROUTINE OVERWRITES * IGCD ON RETURN AND DESTROYS THE CONTENTS OF IX. * * GLOBAL VARIABLES: * * 1. IX - N INTEGER ARRAY - COEFFICIENTS OF IX * 2. IGCD - N INTEGER ARRAY - OUTPUT IGCD COEFFICIENTS * 3. N - INTEGER - NUMBER OF ARRAY ELEMENTS * 4. IBASE - INTEGER BASE OF COEFFICIENTS * 5. IW1,IW2,JU,JV - MP WORK SPACE * ******************************************************************** * INTEGER IX(N),IGCD(N),IW1(N),IW2(N),JU(N),JV(N),IBASE *------------------------------------------------------------------- * FIND LARGER OF IX AND IY AND TAKE THEIR ABS VALUE *------------------------------------------------------------------- DO 10, I = 1, N IW1(I) = 0 IW2(I) = 0 JU(I) = 0 JV(I) = 0 10 CONTINUE CALL MPMAX (IX,IGCD,N,IN) IF ( IN .LE. 1 ) THEN CALL MPCOPY (IW1,IX,N) CALL MPCOPY (IW2,IGCD,N) ELSE CALL MPCOPY (IW1,IGCD,N) CALL MPCOPY (IW2,IX,N) ENDIF IW1(1) = IABS(IW1(1)) IW2(1) = IABS(IW2(1)) *------------------------------------------------------------------- * EUCLIDEAN ALGORITHM *------------------------------------------------------------------- 20 CONTINUE IF ( IW2(1).EQ.1 .OR. IW2(1).EQ.0 ) THEN GO TO 30 ELSEIF ( IW2(2).EQ.1 .AND. IW2(1).EQ.2 ) THEN IW1(1) = 2 IW1(2) = 1 GO TO 30 ELSE CALL MPDIV (IW1,IW2,IX,IGCD,N,IBASE,JU,JV) CALL MPCOPY (IW1,IW2,N) CALL MPCOPY (IW2,IGCD,N) ENDIF GO TO 20 *------------------------------------------------------------------- * SWAP IGCD WITH IW1 ONCE IGCD IS FOUND *------------------------------------------------------------------- 30 CONTINUE CALL MPCOPY (IGCD,IW1,N) RETURN END * ******************************************************************** * SUBROUTINE MPLCM (IX,LCM,N,IBASE,MULT,IW1,IW2,JU,JV) * * FINDS THE LEAST COMMON MULTIPLE OF MP INTEGERS IX,LCM * FROM THE EUCLIDEAN ALGORITHM AND * IX*LCM = GCD(IX,LCM)*LCM(IX,LCM) * THE ROUTINE OVERWRITES LCM ON RETURN AND DESTROYS * THE CONTENTS OF IX. * * GLOBAL VARIABLES: * 1. IX - N INTEGER ARRAY - COEFFICIENTS OF IX * 2. IY - N INTEGER ARRAY - COEFFICIENTS OF IY * 3. LCM - N INTEGER ARRAY - OUTPUT LCM COEFFICIENTS * 4. N - INTEGER - NUMBER OF ARRAY ELEMENTS * 5. IBASE - INTEGER - NUMBER BASE OF IX,IY,IZ * ******************************************************************** * INTEGER IX(N),LCM(N),MULT(N),IW1(N),IW2(N),JU(0:N),JV(0:N), * IBASE * CALL MPMULT (IX,LCM,MULT,N,IBASE) CALL MPGCD (IX,LCM,N,IBASE,IW1,IW2,JU,JV) CALL MPDIV (MULT,LCM,IW1,IW2,N,IBASE,JU,JV) CALL MPCOPY (LCM,IW1,N) RETURN END * ******************************************************************** * SUBROUTINE MPMAKE (IN,IOUT,N,JBASE) * * TO CONVERT ON WORD IN TO MULTIPLE PRECISION OUTPUT WITH * INTEGER RADIX JBASE. THE RESULT IS STORED IN THE * INTEGER ARRAY IOUT. ASSUMES JBASE**2 IS WITHIN THE * WORD SIZE BUT THAT JBASE**3 IS NOT. * * GLOBAL VARIABLES: * * 1. IN - INTEGER - INPUT INTGER * 2. IOUT - N INTEGER ARRAY - OUTPUT EQUAL TO IN * 3. N - INTEGER - SIZE OF IOUT * 4. JBASE - RADIX OF IOUT * ******************************************************************** * INTEGER IOUT(N),IN,JBASE * IOUT(1) = 4 IF ( IN .LT. 0 ) IOUT(1) = -4 INN = IABS(IN) IOUT(4) = MOD(INN,JBASE) INN = (INN-IOUT(4)) / JBASE IOUT(3) = MOD(INN,JBASE) IOUT(2) = (INN-IOUT(3)) / JBASE IF (IOUT(2) .EQ. 0) CALL MPPACK (IOUT,N) RETURN END * ******************************************************************** * SUBROUTINE MPMAX (IX,IY,N,IN) * * TO FIND MAXIMUM OF ABS(IX),ABS(IY) * * GLOBAL VARIABLES: * * 1. IX - N INTEGER ARRAY - COEFFIENTS OF IX * 2. IY - N INTEGER ARRAY - COEFFIENTS OF IY * 3. N - INTEGER - NUMBER OF ARRAY ELEMENTS * 4. IN - INTEGER - OUTPUT AS * IF IX = IY THEN IN = 0 * IF IX > IY THEN IN = 1 * IF IY > IX THEN IN = 2 * ******************************************************************** * INTEGER IX(N),IY(N),IN * IN = 0 NX = IABS(IX(1)) NY = IABS(IY(1)) IF ( NX .GT. NY ) THEN IN = 1 GO TO 20 ELSEIF ( NY .GT. NX ) THEN IN = 2 GO TO 20 ELSE I = 2 IN = 0 10 IF (IX(I) .EQ. IY(I)) THEN I = I + 1 ELSEIF (IX(I) .GT. IY(I)) THEN IN = 1 GO TO 20 ELSE IN = 2 GO TO 20 ENDIF IF (I .LE. NX) GO TO 10 ENDIF 20 RETURN END * ******************************************************************** * SUBROUTINE MPMULT (IX,IY,IZ,N,IBASE) * * TO MULTIPLY MP INTEGERS IX * IY = IZ. * * GLOBAL VARIABLES: * * 1. IX - N INTEGER ARRAY - COEFFIENTS OF IX * 2. IY - N INTEGER ARRAY - COEFFIENTS OF IY * 3. IZ - N INTEGER ARRAY - COEFFIENTS OF IZ * 4. N - INTEGER - NUMBER OF ARRAY ELEMENTS * 5. IBASE - INTEGER - NUMBER BASE OF IX,IY,IZ * * REF: D.E. KNUTH, 'the art of computer programming', * SEC. 4.3.1 * ******************************************************************** * INTEGER IX(N),IY(N),IZ(N),IBASE *-------------------------------------------------------------------- * INITIALIZE *-------------------------------------------------------------------- NX = IABS(IX(1)) NY = IABS(IY(1)) IF ( NX .EQ. 0 .OR. NY .EQ. 0 ) THEN IZ(1) = 0 RETURN ENDIF NZ = NX + NY DO 10, I = 1, NZ IZ(I) = 0 10 CONTINUE IZ(1) = NZ-1 IF ( IX(1)*IY(1) .LT. 0 ) IZ(1) = -IZ(1) *-------------------------------------------------------------------- * MULTIPLY *-------------------------------------------------------------------- DO 30, J = NY-1, 1, -1 K = 0 DO 20, I = NX-1, 1, -1 INT = IX(I+1)*IY(J+1)+IZ(I+J+1)+K IZ(I+J+1) = MOD(INT,IBASE) K = (INT-IZ(I+J+1))/IBASE 20 CONTINUE IZ(J+1) = K 30 CONTINUE IF (IZ(2) .EQ. 0) CALL MPPACK(IZ,N) RETURN END * ******************************************************************** * SUBROUTINE MPOUT (IX,N,IUNIT,KFORM) * * OUTPUTS MP INTEGERS TO IUNIT. RETURNS IX UNCHANGED. * * GLOBAL VARIABLES: * * 1. IX - N INTEGER ARRAY - MP INTEGER TO OUTPUT * 2. N - INTEGER - NUMBER OF ELEMENTS IN IX * 3. IUNIT - INTEGER - UNIT NUMBER OF OUTPUT FILE * 4. KFORM - FLAG INDICATING PRINT MODE WHERE: * = (20*SPACE) + DIGITS * AND WHERE * SPACE = 0 - NO SPACES BETWEEN WORDS * = 1 - ONE SPACE BETWEEN WORDS * DIGITS = LOG10(BASE) WHERE THE CURRENTLY * SUPPORTED VALUES ARE 3-10 * ******************************************************************** * INTEGER IX(N) * ISIGN = 1 IF (IX(1) .LT. 0) ISIGN = -1 IX(2) = IX(2)*ISIGN * IF ( KFORM .LE. 11 ) THEN IF ( KFORM .EQ. 5 ) THEN WRITE (IUNIT,1005) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 6 ) THEN WRITE (IUNIT,1006) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 7 ) THEN WRITE (IUNIT,1007) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 4 ) THEN WRITE (IUNIT,1004) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 8 ) THEN WRITE (IUNIT,1008) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 9 ) THEN WRITE (IUNIT,1009) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 3 ) THEN WRITE (IUNIT,1003) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 10 ) THEN WRITE (IUNIT,1010) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 11 ) THEN WRITE (IUNIT,1011) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 1 ) THEN WRITE (IUNIT,1001) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 2 ) THEN WRITE (IUNIT,1002) (IX(I), I = 2, IABS(IX(1))) ENDIF ELSE IF ( KFORM .EQ. 25 ) THEN WRITE (IUNIT,1025) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 26 ) THEN WRITE (IUNIT,1026) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 27 ) THEN WRITE (IUNIT,1027) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 24 ) THEN WRITE (IUNIT,1024) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 28 ) THEN WRITE (IUNIT,1028) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 29 ) THEN WRITE (IUNIT,1029) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 23 ) THEN WRITE (IUNIT,1023) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 30 ) THEN WRITE (IUNIT,1030) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 31 ) THEN WRITE (IUNIT,1031) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 21 ) THEN WRITE (IUNIT,1021) (IX(I), I = 2, IABS(IX(1))) ELSEIF ( KFORM .EQ. 22 ) THEN WRITE (IUNIT,1022) (IX(I), I = 2, IABS(IX(1))) ENDIF ENDIF * IX(2) = IX(2)*ISIGN RETURN *------------------------------------------------------------------- * SPACE = 0 FORMAT STATEMENTS, NO SPACE *------------------------------------------------------------------- 1001 FORMAT (5X, I2, 65I1, /, 65 (7X,65I1,/)) 1002 FORMAT (5X, I3, 32I2.2, /, 62 (7X,30I2.2,/)) 1003 FORMAT (5X, I5, 15I3.3, /, 62 (7X,16I3.3,/)) 1004 FORMAT (5X, I6, 12I4.4, /, 77 (7X,13I4.4,/)) 1005 FORMAT (5X, I7, 10I5.5, /, 91 (7X,11I5.5,/)) 1006 FORMAT (5X, I8, 8I6.6, /, 112 (7X,9I6.6,/)) 1007 FORMAT (5X, I9, 7I7.7, /, 125 (7X,8I7.7,/)) 1008 FORMAT (5X, I10, 6I8.8, /, 143 (7X,7I8.8,/)) 1009 FORMAT (5X, I11, 5I9.9, /, 167 (7X,6I9.9,/)) 1010 FORMAT (5X, I12, 5I10.10, /, 167 (7X,6I10.10,/)) 1011 FORMAT (5X, I13, 5I11.11, /, 167 (7X,6I11.11,/)) *------------------------------------------------------------------- * SPACE = 1 FORMAT STATEMENTS, ONE SPACE *------------------------------------------------------------------- 1021 FORMAT (5X, I3, 31I2.1, /, 55 (7X,32I2.1,/)) 1022 FORMAT (5X, I4, 19I3.2, /, 62 (7X,20I3.2,/)) 1023 FORMAT (5X, I5, 15I4.3, /, 62 (7X,14I4.3,/)) 1024 FORMAT (5X, I6, 12I5.4, /, 77 (7X,13I5.4,/)) 1025 FORMAT (5X, I7, 10I6.5, /, 91 (7X,11I6.5,/)) 1026 FORMAT (5X, I8, 8I7.6, /, 112 (7X,9I7.6,/)) 1027 FORMAT (5X, I9, 7I8.7, /, 125 (7X,8I8.7,/)) 1028 FORMAT (5X, I10, 6I9.8, /, 143 (7X,7I9.8,/)) 1029 FORMAT (5X, I11, 5I10.9, /, 167 (7X,6I10.9,/)) 1030 FORMAT (5X, I12, 5I11.10, /, 167 (7X,6I11.10,/)) 1031 FORMAT (5X, I13, 5I12.11, /, 167 (7X,6I12.11,/)) END * ******************************************************************** * SUBROUTINE MPPACK (IX,N) * * TO STRIP LEADING TERMS OF IX THAT ARE ZERO, I.E. RETURNS * IX BUT WITH FEWER TERMS * * GLOBAL VARIABLES: * * 1. IX - N INTEGER ARRAY - COEFFIENTS OF IX * 2. N - INTEGER - NUMBER OF ARRAY ELEMENTS * ******************************************************************** * INTEGER IX(N) * NX = IABS(IX(1)) ISIGN = 1 IF (IX(1) .LT. 0) ISIGN = -1 * 10 CONTINUE IF( IX(2) .EQ. 0 ) THEN DO 20, J = 2, NX-1, 1 IX(J) = IX(J+1) 20 CONTINUE IX(NX) = 0 NX = NX-1 ELSE GO TO 30 ENDIF IF (NX .GT. 1) GO TO 10 30 IX(1) = NX*ISIGN IF (NX .EQ. 1) IX(1) = 0 RETURN END * ******************************************************************** * SUBROUTINE MPPWR (IX,NPWR,IY,N,IBASE,IW1,IW2,IW3) * * EXPONENTIATES MP INTEGER IX TO NPWR (IX**NPWR). * GIVES A MP INTEGER RESULT IN ARRAY IY. USES THE * RIGHT TO LEFT BINARY ALGORITHM. * * GLOBAL VARIABLES: * * 1. IX - N INTEGER ARRAY - COEFFIENTS OF IX * 2. NPWR - POWER OF IX DESIRED, NPWR > 0 * 3. IY - N INTEGER ARRAY - COEFFIENTS OF IY * 4. N - NUMBER OF ARRAY ELEMENTS * 5. IBASE - NUMBER BASE OF IX,IY,IZ * 6. IW1-3 - N INTEGER ARRAYS - WORK SPACE * * REF: D.E. KNUTH, 'the art of computer programming', * SEC. 4.6.3 * ******************************************************************** * INTEGER IX(N),IY(N),BASE,NPWR,IW1(N),IW2(N),IW3(N) *-------------------------------------------------------------------- * INTIALIZE *-------------------------------------------------------------------- IF ( NPWR .LE. 0 ) RETURN IY(1) = 2 IY(2) = 1 NN = NPWR CALL MPCOPY (IW1,IX,N) *-------------------------------------------------------------------- * HALVE N AND DETERMINE WHETHER IT WAS EVEN OR ODD. IF ODD * THEN MULTIPLY IY BY IW1. SQUARE IW1. REPEAT UNTIL * NN = 0. *-------------------------------------------------------------------- 10 CONTINUE KK = MOD(NN,2) NN = (NN - KK) / 2 IF ( KK .EQ. 1 ) THEN CALL MPMULT (IW1,IY,IW2,N,BASE) CALL MPCOPY (IY,IW2,N) ENDIF IF ( NN .EQ. 0 ) GO TO 20 CALL MPCOPY (IW2,IW1,N) CALL MPMULT (IW1,IW2,IW3,N,BASE) CALL MPCOPY (IW1,IW3,N) GO TO 10 20 CONTINUE RETURN END * ******************************************************************** * SUBROUTINE MPREAD (IN,N,IUNIT,JBASE,IW1) * * TO READ MULTIPLE PRECISION INPUT AND CONVERT IT TO INTEGER * RADIX JBASE. THE RESULT IS STORED IN THE INTEGER ARRAY IN. * * THE PROCEDURE TO READ THE DATA IS TO READ EACH INTEGER * AS A CHARACTER AND CONVERT ONE DIGIT AT A TIME TO * INTEGER FORMAT. * * WARNING: * * THIS ROUTINE RELIES ON THE FORTRAN COMPILER FUNCTION ICHAR * TO INCREMENT APLHA-NUMERICALLY I.E. IF ICHAR('0') = J THEN * TO EXECUTE PROPERLY IT REQUIRES ICHAR('1') = J+1, * ICHAR('2') = J+2, ... , ICHAR('9') = J+9. * * GLOBAL VARIABLES: * * 1. IN - N INTEGER ARRAY - OUTPUT MP INTEGER RADIX JBASE * 2. N - PHYSICAL DIMENSION OF IN * 3. IUNIT - INPUT FILE UNIT NUMBER. * 4. JBASE - RADIX OF IN ON OUTPUT * 5. IW1 - N INTEGER ARRAY - WORK SPACE IN MPBASE * ******************************************************************** * PARAMETER ( MAXC = 1000) CHARACTER*1 A(80) INTEGER II(MAXC),IN(N),IW1(N) *-------------------------------------------------------------------- * INTIALIZE - GET INPUT AS CHARACTERS *-------------------------------------------------------------------- DO 10, I = 1, N IN(I) = 0 10 CONTINUE I = 1 ISIGN = 1 INDEX = 0 IBASE = 4 IF ( MOD(JBASE,10) .EQ. 0 ) IBASE = INT(ALOG10(FLOAT(JBASE+1))) ISHFT = ICHAR('0') *------------------------------------------------------------------- * CONVERT CHARACTERS TO INTEGERS ONE AT A TIME *------------------------------------------------------------------- 12 CONTINUE READ (IUNIT,1000) (A(KKK),KKK=1,80) DO 20 K = 1,80 IAK = ICHAR(A(K)) - ISHFT IF ( A(K) .EQ. '$') THEN I = I - 1 GO TO 30 ELSEIF ( A(K) .EQ. '-') THEN INDEX = 1 ISIGN = -1 I = I + 1 ELSEIF ( IAK.GE.0 .AND. IAK.LE.9 ) THEN II(I) = ICHAR(A(K)) - ISHFT I = I + 1 ENDIF 20 CONTINUE GO TO 12 *------------------------------------------------------------------- * CONVERT SINGLE INTEGERS TO BASE 10**(IBASE) NUMBERS *------------------------------------------------------------------- 30 CONTINUE IF ( ISIGN .EQ. -1 ) I = I-1 JJ = MOD(I,IBASE) IF ( JJ.EQ.0 ) THEN KOUNT = 0 JJ = I/IBASE ELSE DO 40, K = 1, JJ INDEX = INDEX+1 IN(2) = IN(2) + (II(INDEX)*(10**(JJ-K))) 40 CONTINUE KOUNT = 1 JJ = (I-JJ)/IBASE ENDIF DO 60, K = 1, JJ, 1 KOUNT = KOUNT + 1 DO 50, L = 0, IBASE-1 INDEX = INDEX + 1 IN(KOUNT+1) = IN(KOUNT+1) + (II(INDEX)*(10**(IBASE-1-L))) 50 CONTINUE 60 CONTINUE IN(1) = ISIGN * (KOUNT + 1) *------------------------------------------------------------------- * CONVERT RADIX 10**(IBASE) NUMBERS TO RADIX JBASE NUMBERS *------------------------------------------------------------------- NBASE = 10**IBASE IF ( JBASE .NE. NBASE ) CALL MPBASE (IN,N,NBASE,JBASE,IW1) IF ( IN(2) .EQ. 0 ) CALL MPPACK (IN,N) RETURN 1000 FORMAT(80A1) END * ******************************************************************** * SUBROUTINE MPSUB (IX,IY,IZ,N,BASE) * * TO SUBTRACT MP NON NEGATIVE INTEGERS IX - IY = IZ, WHERE * IX > IY * * GLOBAL VARIABLES: * * 1. IX - N INTEGER ARRAY - COEFFIENTS OF IX * 2. IY - N INTEGER ARRAY - COEFFIENTS OF IY * 3. IZ - N INTEGER ARRAY - COEFFIENTS OF IZ * 4. N - INTEGER - NUMBER OF ARRAY ELEMENTS * 5. BASE - INTEGER - NUMBER BASE OF IX,IY,IZ * * REF: D.E. KNUTH, 'the art of computer programming', * SEC. 4.3.1 * ******************************************************************** * INTEGER IX(N),IY(N),IZ(N),BASE * NX = IABS(IX(1)) NY = IABS(IY(1)) NZ = MAX0(NX,NY) MZ = MIN0(NX,NY) K = 0 DO 10, J = 0, MZ-2 INT = IX(NX-J) - IY(NY-J) + K IF (INT .LT. 0) THEN IZ(NZ-J) = INT + BASE K = -1 ELSE IZ(NZ-J) = INT K = 0 ENDIF 10 CONTINUE IF ( NZ .NE. MZ ) THEN IF ( NZ .EQ. NX ) THEN IZ(NZ-MZ+1) = IX(NZ-MZ+1) + K DO 20, J= NZ-MZ, 2, -1 IZ(J) = IX(J) 20 CONTINUE ELSE IZ(NZ-MZ+1) = IY(NZ-MZ+1) + K DO 30, J= NZ-MZ, 2, -1 IZ(J) = IY(J) 30 CONTINUE ENDIF ENDIF IZ(1) = NZ IF (IZ(2) .EQ. 0) CALL MPPACK(IZ,N) RETURN END * ********************************************************************* * SUBROUTINE PRMS(NPRMS,IPRMS,NPMAX,IWORK,LNGTH) * * THIS SUBROUTINE COMPUTES THE NPRMS PRIMES SMALLER THAN * NPMAX AND STORES THEM IN THE ARRAY IWORK IN DESCENDING * ORDER. * * LAST CHANGE: 8/27/88 * * THE PARAMETERS OF PRMS HAVE THE FOLLOWING MEANING: * * NPRMS: ON INPUT, THE NUMBER OF PRIMES THAT ARE TO BE FOUND. * ON OUTPUT, THE NUMBER OF PRIMES THAT HAVE BEEN FOUND. * THE ONLY RASON WHY FEWER PRIMES THAN REQUESTED MIGHT BE * FOUND IS THAT THERE ARE NOT SUFFICIENTLY MANY PRIMES * SMALLER THAN NPMAX. * * IPRMS: A ONE DIMENSIONAL INTEGER ARRAY OF SUFFICIENT SIZE THAT * ON OUTPUT CONTAINS THE REQUIRED PRIMES IN DESCENDING * ORDER * * NPMAX: THE MAXIMUM SIZE OF THE PRIMES TO BE FOUND. NPMAX IS NOT * CHANGED BY THE SUBROUTINE AND MAY BE EVEN OR ODD. * * IWORK: AN INTEGER WORK ARRAY THAT IS TO BE AS LARGE AS POSSIBLE. * THE LARGER IT IS THE MORE EFFICIENT WILL BE THE GENERATION * OF THE PRIMES. THIS SUBROUTINE IS BASED ON AN INTEGER * WORK ARRAY - RATHER THAN A LOGICAL ONE - SINCE ITS ANTICIPATED * USE IS IN AN EVIRONMENT WHERE LARGE INTEGER WORK SPACE IS AVAI * BUT NO ADDITIONAL MEMORY CAN BE SUPPLIED. * * LNGTH: THE PHYSICAL LENGTH OF IWORK. * * ERROR CONDITIONS: * * THE SUBROUTINE DOES NOT USE ANY IO. THE PARAMETER NPRMS WILL EQUAL * ZERO ON OUTPUT IF ONE OF THE FOLLOWING CONDITIONS OCCURS: * * 1. NPMAX < 2 ON INPUT. * * 2. LNGTH < 2 ON INPUT * * THE PRIME GENERATION IS BASED ON THE SIEVE OF ERATOSTHENES. * ********************************************************************* * INTEGER NPRMS,IPRMS(1),NPMAX,IWORK(1),LNGTH * * CHECK FOR INCONSISTENCIES IN INPUT DATA: * IF (NPMAX.LT.2.OR.LNGTH.LT.2) THEN NPRMS = 0 RETURN ENDIF * * MAKE SURE UPPER BOUND IS ODD * NPM = NPMAX - 1 + MOD(NPMAX,2) * * SET UP CALL FOR SIEVING ROUTINE. ORDINARILY, IT WILL BE CALLED * JUST ONCE * NP = NPRMS NFOUND = 1 NCNT = 0 * * CALL SIEVING ROUTINE * 100 CONTINUE CALL PSIEVE(NP,IPRMS(NFOUND),NPM,IWORK,LNGTH) NCNT = NCNT + NP IF (NP.EQ.0.OR.NCNT.EQ.NPRMS) GO TO 200 NPM = IPRMS(NCNT) - 2 NP = NPRMS - NCNT NFOUND = NCNT + 1 GO TO 100 200 CONTINUE IF (NCNT.LT.NPRMS.AND.(IPRMS(NCNT).EQ.3.OR.NCNT.EQ.0)) THEN NCNT = NCNT + 1 IPRMS(NCNT) = 2 ENDIF IF (NPRMS.NE.NCNT) NPRMS = NCNT END * SUBROUTINE PSIEVE(NPRMS,IPRMS,NPMAX,IWORK,LNGTH) * * CALLING LIST IS IDENTICAL TO CALLING LIST FOR PRMS * INTEGER NPRMS,IPRMS(1),NPMAX,IWORK(1),LNGTH * * THE FIRST FOUR PRIMES ARE STORED IN NFIRST(4) INTEGER NFIRST(4) NFIRST(1) = 3 NFIRST(2) = 5 NFIRST(3) = 7 NFIRST(4) = 11 * * SET A SUITABLE SEARCH INTERVAL * USE THE FACT THAT THE DENSITY OF PRIMES IS APPROXIMATELY EQUAL * TO 1/LN(NPMAX), AND USE 1.1 TIMES THE INTERVAL FOR SAFETY. * IROOT = IFIX(SQRT(FLOAT(NPMAX))) IROOT = IROOT - 1 + MOD(IROOT,2) NSRCH = MAX(100,IFIX(FLOAT(NPRMS)*ALOG(FLOAT(NPMAX))*1.1)) NSRCH = MIN(NSRCH,NPMAX-2,2*LNGTH,NPMAX-IROOT) NSRCH = NSRCH + MOD(NSRCH,2) NSRCH2 = NSRCH/2 * * IWORK(I) WILL BE 1 IF NPRMS - 2*(I-1) IS PRIME * AND 0 OTHERWISE * * FIRST PRESUME ALL NUMBERS IN THE RANGE ARE PRIME * DO 100 I = 1,NSRCH2 IWORK(I) = 1 100 CONTINUE * * LABEL ALL LOCATIONS 0 THAT ARE MULTIPLES OF 3,5,7, OR 11. * DO 300 K = 1,4 IINC = NFIRST(K) ISTART = NPMAX - MOD(NPMAX,IINC) IF (MOD(ISTART,2).EQ.0) ISTART = ISTART - IINC ISTART = (NPMAX - ISTART + 2)/2 ILAST = MIN(NSRCH2,(NPMAX-3*IINC)/2+1) DO 200 I = ISTART,ILAST,IINC IWORK(I) = 0 200 CONTINUE 300 CONTINUE * * NOW GO THROUGH ODD NUMBER GREATER THAN 11 AND IGNORE MULTIPLES OF * 3,5,7, OR 11. * DO 500 IINC = 13,IROOT,2 IF (MOD(IINC,3).EQ.0.OR.MOD(IINC,5).EQ.0.OR.MOD(IINC,7) X .EQ.0.OR.MOD(IINC,11).EQ.0) GO TO 500 ISTART = NPMAX - MOD(NPMAX,IINC) IF (MOD(ISTART,2).EQ.0) ISTART = ISTART - IINC ISTART = (NPMAX - ISTART + 2)/2 ILAST = MIN(NSRCH2,(NPMAX-3*IINC)/2+1) DO 400 I = ISTART,ILAST,IINC IWORK(I) = 0 400 CONTINUE 500 CONTINUE * * WE HAVE FOUND THE PRIMES, NOW STORE THEM * NPFND = 0 DO 600 I = 1,NSRCH2 IF (NPFND.EQ.NPRMS) GO TO 700 IF (IWORK(I).EQ.1) THEN NPFND = NPFND + 1 IPRMS(NPFND) = NPMAX-2*(I-1) ENDIF 600 CONTINUE 700 CONTINUE NPRMS = NPFND RETURN END * ********************************************************************* * SUBROUTINE RADIXF(IC,N,IA,IWRK,NWRK) * * GIVEN A MIXED-RADIX INTEGER COMPUTED WITH RADIXM THIS ROUTINE * COMPUTES THE COEFFICIENTS OF ITS FIXED-RADIX REPRESENTATION. * BOTH FORMS ARE STORED IN THE IWRK ARRAY. THE INPUT MIXED * RADIX FORM BEGINS AT LOCATION IC, AND THE OUTPUT FIXED RADIX * FORM BEGINS AT LOCATION IA. THE OUTPUT IS IN MP FORM. * * GLOBAL VARIABLES: * * 1. IC - MIXED RADIX COEFFICIENTS ADDRESS * 2. N - INTEGER - # OF COEFS. OF MIXED RADIX FORM * 3. IA - ADDRESS OF FIXED RADIX COEFICIENTS * 4. IWRK - NWRK INTEGER ARRAY - SEE DRIVER * * COMMON VARIABLES USED: * 1. IAPRM - ACTUAL PRIME ADDRESS * 2. IBASE - OUTPUT BASE * ********************************************************************* * INTEGER IWRK(NWRK),IC,IA,N INTEGER IQ,IQTMP,IPP,L,ISIGN COMMON /WPARM/ MINFO,IHPRM,IROW,ICOL,IDEP,IPIV,LST,IPRM,IPIN, * IAPRM,KODE,IRHS,JWRK,KNUM,KGCD,IMAT,KUN13,IBASE, * ITER,ICPRM,IDIAG,ITRM,IBASIC,ITRI,INULL,IRANK,IZERO,LMAX, * NPRME,MHZ,ISOLH,KUN12 L = 1 IWRK(IA+1) = 0 *------------------------------------------------------------------- * AT THIS STAGE IWRK(IA+1) - IWRK(IA+L) IS THE FIXED-RADIX REP. * OF IWRK(IC+N-I+2) - IWRK(IC+N). COMPUTE THE FIRST L COEFS. * OF THE FIXED-RADIX REP. INCLUDING THE NEXT MIXED RADIX COEF. *------------------------------------------------------------------- DO 30, I = 1, N IPP = IWRK(IAPRM+N-I+1) IQ = IWRK(IC+N-I+1) DO 10, J = 1, L IQTMP = IWRK(IA+J)*IPP + IQ IQ = IQTMP/IBASE IWRK(IA+J) = IQTMP - IQ*IBASE 10 CONTINUE *------------------------------------------------------------------- * REPEAT ... CONVERT NEXT TERM TO FIXED-RADIX FORM. *------------------------------------------------------------------- 20 CONTINUE IF ( IQ .NE. 0 ) THEN L = L + 1 IQTMP = IQ/IBASE IWRK(IA+L) = IQ - IQTMP*IBASE IQ = IQTMP GO TO 20 ENDIF 30 CONTINUE *------------------------------------------------------------------- * REORDER THE COEFFICIENTS. *------------------------------------------------------------------- IF (L.GE.2) THEN DO 40, I = 1, L/2 IQTMP = IWRK(IA+I) IWRK(IA+I) = IWRK(IA+L-I+1) IWRK(IA+L-I+1) = IQTMP 40 CONTINUE ENDIF *------------------------------------------------------------------- * REPRESENT IN MP FORM *------------------------------------------------------------------- ISIGN = 1 IF ( IWRK(IA+1) .LT. 0 ) ISIGN = -1 DO 50, I = L, 1, -1 IWRK(IA+I+1) = IABS(IWRK(IA+I)) 50 CONTINUE IWRK(IA+1) = (L+1)*ISIGN IF ( L .EQ. 0 ) IWRK(IA+1) = 0 RETURN END * ******************************************************************** * SUBROUTINE RADIXM (ILOC,IRES,TEST,IWRK,NWRK) * * GIVEN THE RESIDUUM IRES MODULO ICPRM AND THE FIRST ITER-1 * COEFFICIENTS IN THE MIXED-RADIX REPRESENTATION OF AN * INTEGER, THIS ROUTINE COMPUTES THE ITER-TH COEFFICIENT. * * GLOBAL VARIABLES: * * 1. ILOC - ADDRESS OF PREVIOUS COEFFICIENTS IN IWRK * 2. IRES - RESIDUE OF NEW COEFFICIENT * 3. TEST - LOGICAL - NOT CHANGED IF NEW RADIX IS ZERO * ELSE TEST IS SET .FALSE.. * 4. IWRK - NWRK INTEGER ARRAY - SEE DRIVER * * COMMON VARIABLES USED: * * 1. IAPRM - ACTUAL PRIME LOCATIONS * 2. ICPRM - CURRENT PRIME * 3. IPIN - PRIME INVERSES * 4. ITER - ITERATION COUNTER * ******************************************************************** * INTEGER IWRK(NWRK),ILOC,IRES INTEGER NTEMP,I1 LOGICAL TEST COMMON /WPARM/ MINFO,IHPRM,IROW,ICOL,IDEP,IPIV,LST,IPRM,IPIN, * IAPRM,KODE,IRHS,JWRK,KNUM,KGCD,IMAT,KUN13,IBASE, * ITER,ICPRM,IDIAG,ITRM,IBASIC,ITRI,INULL,IRANK,IZERO,LMAX, * NPRME,MHZ,ISOLH,KUN12 * NTEMP = IRES IF ( ITER .GT. 1 ) THEN NTEMP = IWRK(ILOC+ITER-1) DO 10 I = 1, ITER-2, 1 I1 = ITER - 1 - I NTEMP = MOD(NTEMP*IWRK(IAPRM+I1)+IWRK(ILOC+I1),ICPRM) 10 CONTINUE NTEMP = MOD(IWRK(IPIN+ITER)*MOD(IRES-NTEMP,ICPRM),ICPRM) ENDIF IF ( NTEMP .NE. 0 ) TEST = .FALSE. I1 = (ICPRM+1)/2 IWRK(ILOC+ITER) = NTEMP - (NTEMP/I1)*ICPRM RETURN END * ******************************************************************** * SUBROUTINE RSDANA (IPARM,IHTBL,NSTOR,IWRK,NWRK,HWRK,NH) * * DRIVER ROUTINE * * THIS ROUTINE ANALYZES LARGE SPARSE INTEGER SYSTEMS OF LINEAR * EQUATIONS USING RESIDUE TECHNIQUES. IT REQUIRES AN INPUT * MATRIX FILE DESCRIBED BELOW AND REMAINS IN-CORE EXECPT * TO READ THE MATRIX FILE. * * EXPLICITLY, BASED ON INPUT PARAMETERS, IT FINDS A PARTICULAR * SOLUTION THAT IS A LEAST SQUARES SOLUTION, ALTHOUGH NOT THE * MINIMUM NORM LEAST SQUARES SOLUTION, A BASIS FOR THE NULL * SPACE AND, THE RANK AND ROW DEPENDENCIES OF THE EQUATION * * AX = B * WHERE: * A IS AN NXM INTEGER MATRIX * X IS THE MXP UNKNOWN MATRIX * B IS AN NXP INTEGER RIGHT HAND SIDE (RHS) * * THIS ROUTINE ASSUMES NO PARTICULAR STRUCTURE ON MATRIX A * OR THAT A IS IN FACT SPARSE, HOWEVER, IT HAS HIGH OVERHEAD * FOR NON-SPARSE PROBLEMS AND AS SUCH IS NOT RECOMMENDED. * IT ALLOWS ARBITRARILY LARGE INTEGERS ON INPUT. STOPPING * CRITERION ARE OPTIONAL BUT NOT GUARANTEED TO BE MET, THE * PROGRAM WILL TERMINATE PREMATURELY IF THE HASHING DATA * STRUCTURE BECOMES 100% FULL. * * THE ROUTINE ITSELF IT IS WRITTEN WITHIN THE FORTRAN-77 * STANDARD AND IS MEANT TO BE FULLY PORTABLE. THE ONLY * SYMBOLS USED IN THE ACTUAL CODE ARE: * 0 1 2 3 4 5 6 7 8 9 * ABCDEFGHIJKLMNOPQRSTUVWXYZ * 'blank' = + - * / ( ) , . $ ' : * * THE ROUTINE WILL TRY TO WRITE TO UNIT 3 IF NO OUTPUT * TERMINAL UNIT NUMBER IS SUPPLIED. * * ******************************************************************** * * VARIABLE LISTS ARE ALL INTEGERS UNLESS OTHERWISE SPECIFIED * ****************************INPUT*********************************** * * * GLOBAL VARIABLES: * * 1. IPARM - 0:15 INTEGER ARRAY - MATRIX/ROUTINE PARAMETERS WHERE * (0) - MAXIMUM NUMBER OF ROWS IN AUGMENTED MATRIX * (1) - MAXIMUM NUMBER OF COLUMNS IN AUGMENTED MATRIX * (2) - MAXIMUM NUMBER OF RIGHT HAND SIDES * (3) - INTEGER INPUT FILE NUMBER * (4) - HARD COPY OUTPUT FILE NUMBER * (5) - INTERMEDIATE OUTPUT FILE NUMBER * (6) - TERMINAL OUTPUT UNIT NUMBER * (7) - INTERMEDIATE OUTPUT FLAG * = 0 THEN NONE * = 1 TERMINAL PROGRESS REPORT * = 2 PIVOT INFO (PASS ONE) * = 3 MATRIX INFO (PASS ONE) * = 4 ALL INFO (CAUTION) * (8) - TASK INDICATOR * = 0 RANK AND ROW DEPENDENCIES ONLY * = 1 RANK AND PARTICULAR SOLUTION * = 2 RANK AND NULL SPACE * = 3 RANK, PARTICULAR SOLUTION AND NULL SPACE * = 4 INTEGER COVNERSION ONLY * (9) - OUTPUT NUMBER BASE * = USUALLY (BUT NOT NECESSARILY) SOME POWER OF * 10 SO THAT ITS SQUARE IS WITHIN COMPUTER WORD * (10) - STOPPING CRITERION ON NUMBER OF PASSES * = -2 MODIFIED HADAMARD BOUND * = -1 ROUTINE DECIDES * = 0 HADAMARD BOUND * = N - SOME POSITIVE INTEGER (USER DECIDES) * (11) - LARGEST INTEGER CAPABLE OF BEING STORED * (12) - MAXIMUM NUMBER OF PRIMES ALLOWED * (13) - OUTPUT STYLE FOR MP INTEGERS * (14) - TYPE OF RATIONAL FILE (1=SP RATIONAL, 2=MP RATIONAL) * TYPE = 0 MEANS INTEGER FILE * (15) - IO NUMBER FOR RATIONAL FILE * * 2. IHTBL - 2 X NSTOR INTEGER ARRAY - HASH STORAGE AREA * * 3. NSTOR - PHYSICAL DIMENSION OF IHTBL AND IT SHOULD BE * AS LARGE AS POSSIBLE. * * 4. IWRK - NWRK INTEGER ARRAY - WORK SPACE * * 5. NWRK - PHYSICAL DIMENSION OF IWRK. THE INTEGER * FUNCTION IWORK IS PROVIDED TO CALCULATE THE REQUIRED * DIMENSION OF IWRK. * * 6. HWRK - NH REAL ARRAY - HADAMARD BOUND WORK SPACE * * 7. NH - PHYSICAL DIMENSION OF HWRK WHERE * NH > MAX(IPARM(0)+IPARM(2),IPARM(1)) * * * INPUT FILE STUCTURE - MPCON OUTPUT STRUCTURE. * NOTE, THE MULTIPLE PRECISION INTEGERS IN THIS FILE ARE * IN A SLIGHTLY DIFFERENT FORMAT THAN IS REQUIRED BY THE * MP ARITHMETIC ROUTINES. TO GET THAT FORMAT SIMPLY * MULTIPLY (AIJ)1 BY NTRM. * * NROW, NCOL, NRHS, NELE, JBASE * I, J, NTRM ] * (AIJ)K , K=1,2,...,NTRM ] - REPEATED NELE TIMES * * WHERE: NROW = # OF ROWS * NCOL = # OF COLUMNS * NRHS = # OF RIGHT HAND SIDES * NELE = # OF NON-ZERO ELEMENTS * I,J = ELEMENT MATRIX INDEX * NTRM = NUMBER OF TERMS RADIX(JBASE) OF AIJ * (AIJ)K = AIJ ELEMENT WHERE K IS THE TERM INDEX * AND (AIJ)1 REPRESENTS THE SIGN OF THE ELEMENT. * THUS AIJ IS * * AIJ = (AIJ)1 * ( (AIJ)2*(JBASE**(NTMR-2)) + * (AIJ)3*(JBASE**(NTMR-3)) + * ... * (AIJ)NTRM-1*JBASE + (AIJ)NTRM ) * * ****************************OUTPUT********************************** * * * 1. IPARM IS OVERWRITTEN ON OUTPUT AS FOLLOWS: * * IPARM: * * (0) - ERROR CONDITION INDICATOR * = 0 THEN NO ERROR * = -1 ERROR ON INPUT * = -2 STOP CRITERION NOT SATISFIED * (1) = LENGTH OF THE HASH TABLE * (2) = HASHING MULTIPLIER * (3) = USEABLE HASH TABLE SIZE * (4) = MATRIX RANK * (5) = COLUMNS OF AUGMENTED MATRIX * (6) = ROWS OF MATRIX + 1 * (7) = MAXIMUM DOUBLE HASHES TRIED * (8) = DIMENSION OF THE NULL SPACE * (9) = TOTAL DOUBLE HASHES TRIED * (10) = TOTAL HASH STORAGES TRIED * (11) = TERMINAL UNIT NUMBER * (12) = ADDRESS OF DEPENDENT ROWS IN IWRK * (13) = ADDRESS OF RHS CONSISTENCY FLAG IN IWRK * * 2. SOLUTIONS AND NULL SPACE BASIS ARE STORED IN THE HASH * TABLE IN MP FORM. * * TO ACCESS THE SOLUTIONS THE USER MUST USE THE HASHING * DATA STRUCTURE. THE INTEGER FUNCTION IRHASH RECOVERS * THE DATA INDICIES IN THE HASH TABLE. AN EXAMPLE * OF ITS USE IS PROVIDED BELOW, AND FOR CONVIENIENCE * IPARM IS OVERWRITTEN TO BE USED IN THE CALL STATEMENT * TO IRHASH. * * THE RIGHT HAND SIDE SOLUTIONS ARE ARE ACCESSED BY * THEIR ROW INDEX AND THEIR AUGMENTED MATRIX COLUMN * INDEX. THE DENOMINATOR FOR ANY RHS IS STORED * IN ((# OF ROWS) + 1) "ROW" LOCATION. FOR EXAMPLE, * IF A IS A NROWxNCOL MATRIX AND HAS 1 RHS, TO ACCESS * THE RESULTS FOR THAT RHS AND STORE THEM IN A 2D ARRAY * CALLED IRHS, THAT HAS THE ROW INDEX AS THE FIRST * DIMENSION AND TERM COUNTER IN THE SECOND DIMENSION, * THE USER WOULD HAVE CODE SIMILAR TO * * DO 20, I = 1, NCOL + 1 * KEY = IRHASH (I,NCOL+1,1,IPARM(1),IHTBL,NSTOR) * IRHS(I,1) = IHTBL(2,KEY) * DO 10, K = 2, IABS(IRHS(I,1)) * KEY = IRHASH (I,NCOL+1,K,IPARM(1),IHTBL,NSTOR) * IRHS(I,K) = IHTBL(2,KEY) * 10 CONTINUE * 20 CONTINUE * * THE NULL SPACE IS ACCESSED THROUGH A COLUMN NUMBER THAT * RANGES OVER 1,2,...,(DIMENSION OF THE NULL SPACE) AND * A ROW INDEX THAT RANGES OVER 1,2,...,(# OF COLUMNS OF A). * THE NULL SPACE BASIS IS SCALED SO THAT ALL ITS ELEMENTS * ARE INTEGERS. THUS NO DENOMINATOR IS REQUIRED. * IT IS ACCESSED WITH SIMILAR CALLS TO IRHASH AS ABOVE. * FOR EXAMPLE, TO GET ALL THE VECTORS IN THE NULL SPACE * FROM THE HASH TABLE AND STORE THEM IN A 3D ARRAY (INULL) * THE USER WOULD NEED CODE SIMILAR TO * * DO 30, J = 1, IPARM(8) * DO 20, I = 1, NCOL * KEY = IRHASH (J,I,1,IPARM(1),IHTBL,NSTOR) * INULL(I,J,1) = IHTBL(2,KEY) * DO 10, K = 2, IABS(INULL(J,I,1)) * KEY = IRHASH (J,I,K,IPARM(1),IHTBL,NSTOR) * INULL(J,I,K) = IHTBL(2,KEY) * 10 CONTINUE * 20 CONTINUE * 30 CONTINUE * * FINALLY, THE MINOR OF THE LEADING SQUARE NON-SINGULAR * MATRIX (I.E. THE DETERMINANT OF AN INVERTIBLE MATRIX) * CAN BE RECOVERED (NCOL+1,1) POSITION AS * * KEY = IRHASH (NCOL+1,1,1,IPARM(1),IHTBL,NSTOR) * IDEN(1) = IHTBL(2,KEY) * DO 10, K = 2, IABS(IDEN(1)) * KEY = IRHASH (NCOL+1,1,K,IPARM(1),IHTBL,NSTOR) * IDEN(J,I,1) = IHTBL(2,KEY) * 10 CONTINUE * ****************************LOCAL*********************************** * * * COMMON /WPARM/ LOCAL VARIABLES: * * IAPRM - PRIME INVERSES ADDRESS IN IWRK * IBASE - OUTPUT NUMBER BASE * IBASIC - PARTICULAR SOLUTION OUTPUT FLAG * ICOL - COLUMN PIVOT VECTOR ADDRESS IN IWRK * ICPRM - CURRENT PRIME * IDEP - DEPENDENT ROW INFORMATION ADDRESS IN IWRK * IDIAG - LEADING MATRIX MINOR RESIDUE * IHPRM - HASH PARAMETER ADDRESS IN IWRK * IMAT - INTERMEDIATE OUTPUT FLAG * INULL - NULL SPACE OUTPUT FLAG * IPIN - PRIME INVERSES ADDRESS IN IWRK * IPIV - ROW/COL NON-ZERO ELEMENT COUNTER ADDRESS IN IWRK * IPRM - PRIMES ADDRESS IN IWRK * IROW - ROW PIVOT VECTOR ADDRESS IN IWRK * ITER - PASS COUNTER * ITRI - TRIANGULAR DECOMPOSITION FLAG * ITRM - TERMINAL UNIT NUMBER * IRANK - RESIDUE RANK FLAG * IRHS - RHS CONSISTENCY FLAG ADDRESS IN IWRK * ISOLH - SOLUTION HASH TABLE FULL FLAG * IZERO - RESIDUE FLAG * JWRK - WORK ADDRESS IN IWRK * KGCD - GREATEST COMMON DIVISOR ADDRESS IN IWRK * KODE - COMMON DENOMINATOR ADDRESS IN IWRK * KNUM - SOLUTION WORK SPACE IN IWRK * KUN13 - INTERMEDIATE OUTPUT FILE NUMBER * KUN12 - OUTPUT FILE NUMBER * LMAX - MAXIMUM NUMBER OF TERMS ON OUTPUT * LST - WORK SPACE ADDRESS IN IWRK * NPRME - MAX NUMBER OF PRIMES * MINFO - MATRIX INFORMATION ADDRESS IN IWRK * MHZ - ZERO POSITION OF THE SOLUTION HASH TABLE * * IWRK WORK SPACE ALLOCATION (INTERNAL): * MINFO * (MINFO+1) = NUMBER OF ROWS IN MATRIX * (MINFO+2) = NUMBER OF COLUMNS IN MATRIX * (MINFO+3) = NUMBER OF RIGHT HAND SIDES * (MINFO+4) = NUMBER OF NONZERO ELEMENTS IN MATRIX FILE * (MINFO+5) = MINIMUM NUMBER OF ROWS OR COLUMNS * (MINFO+6) = MAXIMUM NUMBER OF ROWS OR COLUMNS * (MINFO+7) = RANK APPROXIMATION * (MINFO+8) = ELIMINATION EFFICIENY COUNTER * (MINFO+9) = FILLIN TOTAL * IHPRM * (IHPRM+1) = LENGTH OF THE HASH TABLE * (IHPRM+2) = HASHING MULTIPLIER * (IHPRM+3) = ELIMINATION HASH TABLE SIZE * (IHPRM+4) = SOLUTION HASH TABLE SIZE * (IHPRM+5) = COLUMNS OF AUGMENTED MATRIX * (IHPRM+6) = MAX. ROWS OR (COLUMNS+1) * (IHPRM+7) = ELIMINATION MAX. DOUBLE HASHES TRIED * (IHPRM+8) = SOLUTION MAX. DOUBLE HASHES TRIED * (IHPRM+9) = TOTAL DOUBLE HASHES TRIED * (IHPRM+10) = TOTAL HASH STORAGES TRIED * (IHPRM+11) = TERMINAL UNIT NUMBER * (IHPRM+12) = ORIGINAL VALUE OF NSTOR * (IHPRM+13) = SOLUTION HASH TABLE USE * (IHPRM+14) = NON ZERO ELEMENTS IN RHS * (IHPRM+15) = STOPPING CRITERIA (INTERNAL) * IROW * (IROW) - (IROW + # OF ROWS): ROW PIVOT INDEX * ICOL * (ICOL) - (ICOL + # OF COLS): COLUMN PIVOT INDEX * IDEP * (IDEP) - (IDEP + # OF ROWS): DEPENDENT ROW INDEX * IPIV * (IPIV) - (IPIV + # OF ROWS): FILLIN COUNTER * IRHS * (IRHS) - (IRHS + # OF RHS): RHS CONSISTENCY FLAG * IPRM * (IPRM) - (IPRM + # OF PRIMES): PRIMES STORAGE * IPIN * (IPIN) - (IPIN + # OF PRIMES): INV OF PRIMES * IAPRM * (IAPRM) - (IAPRM + # OF PRIMES): PRIMES USED * KODE * (KODE) - (KODE + LMAX): COMMON DENOM OF SOLUTION * JWRK * (JWRK) - (NWRK): TOTAL WORK SPACE * LST * (LST) - (LST+LMAX): GENERAL WORK SPACE * KNUM * (KNUM) - (KNUM+LMAX): WORK SPACE FOR NUMERATORS * KGCD * (KGCD) - (KGCD+LMAX): WORK SPACE FOR COMMON DIVISOR * IW1-IW2-IW3 * (IW?) - (IW?+LMAX): MP WORK SPACE * * DRIVER LOCAL VARIABLES: * * FILL - REAL - % OF ELIMINATION HASH TABLE FILL * HBASE - REAL - INPUT NUBER BASE * HLOG - REAL - LOG 10 OF HBASE * IDIF - COEF. POWER DIFF. IN PSEUDO-FLOATING MATH * IER - INPUT ERROR CONDITION * IP - 1 INTEGER ARRAY - USED COMPUTING 1 PRIME * ITTOL - TOTAL PASSES REQUIRED * IWTOT - TOTAL WORK SPACE REQUIRED IN IWRK * JBASE - NUMBER BASE OF INPUT FILE * JDEN - ROW OF DENOMINATORS ON OUTPUT * KFORM - FORMAT # OF OUTPUT IN MPOUT * KOUNT - CONVERGENCE COUNTER * LENGTH - AMOUNT OF VOLATILE WORK SPACE * LWRK - VOLATILE WORK SPACE DIVISION * MAXP - PRIME BOUND * NDIS - NUMBER OF DISCARDED PRIMES WHEN RANK DECR. * NPASS - HADAMARD BOUND PASSES REQUIRED * NPASS1 - MODIFIED HADAMARD BOUND PASSES REQUIRED * TEST - LOGICAL FLAG * TOUT - LOGICAL FLAG FOR OUTPUT * ******************************************************************** * COMMON /WPARM/ MINFO,IHPRM,IROW,ICOL,IDEP,IPIV,LST,IPRM,IPIN, * IAPRM,KODE,IRHS,JWRK,KNUM,KGCD,IMAT,KUN13,IBASE, * ITER,ICPRM,IDIAG,ITRM,IBASIC,ITRI,INULL,IRANK,IZERO,LMAX, * NPRME,MHZ,ISOLH,KUN12 * INTEGER IPARM(0:15),IHTBL(2,0:NSTOR),IWRK(NWRK),IP(1) LOGICAL TEST,TOUT REAL HWRK(NH),FILL * * FILL determines percentage of elimnation hash table that will * be filled. A value between 0.3 and 0.7 is recommended. The larger * FILL the faster the progra will run and the more memory will be * required. * PARAMETER ( FILL = 0.5 ) *------------------------------------------------------------------- * INITIALIZE WPARM COMMON BLOCK *------------------------------------------------------------------- IAPRM = 0 IBASE = 0 IBASIC = 0 ICOL = 0 ICPRM = 0 IDEP = 0 IDIAG = 0 IHPRM = 0 IMAT = 0 INULL = 0 IPIN = 0 IPIV = 0 IPRM = 0 IROW = 0 ITER = 0 ITRI = 0 ITRM = 0 IRANK = 0 IRHS = 0 ISOLH = 0 IZERO = 0 JWRK = 0 KGCD = 0 KODE = 0 KNUM = 0 KUN13 = 0 KUN12 = 0 LMAX = 0 LST = 0 NPRME = 0 MINFO = 0 MHZ = 0 *------------------------------------------------------------------- * INPUT PARAMETER CHECKS - I/O UNIT NUMBERS *------------------------------------------------------------------- IER = 0 ITRM = IPARM(6) TEST = .TRUE. IF (ITRM .LE. 0) THEN ITRM = 3 WRITE (ITRM,1180) WRITE (ITRM,1250) IER = 1 ENDIF IF (IPARM(3).LE.0) THEN TEST = .FALSE. ELSE INQUIRE (IPARM(3),OPENED=TEST) ENDIF IF ( .NOT. TEST ) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1100) IER = 1 ENDIF TOUT = .FALSE. IF (IPARM(4) .GT. 0 ) THEN TOUT = .TRUE. INQUIRE (IPARM(4),OPENED=TEST) IF ( .NOT. TEST ) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1110) IER = 1 ENDIF ENDIF *------------------------------------------------------------------- * INPUT PARAMETER CHECKS - MATRIX DIMENSIONS *------------------------------------------------------------------- IF (IPARM(0) .LE. 0) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1220) IER = 1 ENDIF IF (IPARM(1) .LE. 0) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1230) IER = 1 ENDIF IF (IPARM(2).LE.0 .AND. (IPARM(8).EQ.3.OR.IPARM(8).EQ.1)) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1240) IER = 1 ENDIF IPP = IPARM(11) / MAX0(1,(MAX0(IPARM(0),IPARM(1)+1)* * (IPARM(1)+IPARM(2)))) IF ( IPARM(12) .GT. IPP ) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1210) IPP IER = 1 ENDIF IF ( NSTOR .LE. 1 ) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1215) IER = 1 ENDIF *------------------------------------------------------------------- * INPUT PARAMETER CHECKS - OUTPUT PARAMETERS *------------------------------------------------------------------- IPP = 1 IF (IPARM(5).LE.0) THEN TEST = .FALSE. ELSE INQUIRE (IPARM(5),OPENED=TEST) ENDIF IF ( TEST ) IPP = 0 IF (IPARM(7).LT.0 .OR. IPARM(7).GT.4) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1120) IER = 1 ELSEIF (IPARM(7).GT.1 .AND. IPP.EQ.1) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1130) IER = 1 ENDIF IF (IPARM(8).LT.0 .OR. IPARM(8).GT.4) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1140) IER = 1 ENDIF MAXP = INT(SQRT(FLOAT(IPARM(11)))) IF (IPARM(9).GT.MAXP) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1150) MAXP IER = 1 ENDIF IF (IPARM(10) .LT. -2) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1160) IER = 1 ENDIF IF (IPARM(13).LT.0 .OR. IPARM(13).GT.1) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1260) IER = 1 ENDIF *------------------------------------------------------------------- * INPUT PARAMETER CHECKS - WORK ARRAY DIMENSIONS *------------------------------------------------------------------- IF ( NH .LE. MAX(IPARM(0)+IPARM(2),IPARM(1)) ) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1200) MAX(IPARM(0)+IPARM(2),IPARM(1)) IER = 1 ENDIF IPP = IWORK(IPARM) IF ( NWRK .LT. IPP ) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1205) IPP IER = 1 ENDIF *------------------------------------------------------------------- * INPUT PARAMETER CHECKS - RETURN TO MAIN IF ERROR OCCURS *------------------------------------------------------------------- 20 IF ( IER .EQ. 1 ) THEN WRITE (ITRM,1190) RETURN ENDIF *------------------------------------------------------------------- * DATA STRUCTURE INITIALIZATION *------------------------------------------------------------------- 30 CONTINUE DO 40, J = 0, NSTOR IHTBL(1,J) = 0 IHTBL(2,J) = 0 40 CONTINUE DO 50, J = 1, NWRK IWRK(J) = 0 50 CONTINUE *------------------------------------------------------------------- * COMPUTE PARITIONING OF THE IWRK ARRAY *------------------------------------------------------------------- MINFO = 0 IHPRM = 10 IROW = IHPRM + 15 ICOL = IROW + IPARM(0) IDEP = ICOL + IPARM(1) IPIV = IDEP + IPARM(0) IRHS = IPIV + MAX0(IPARM(0),IPARM(1)) IPRM = IRHS + IPARM(2) *------------------------------------------------------------------- * COMPUTE PRIMES - LOAD INTO IWRK *------------------------------------------------------------------- NPRME = IPARM(12) CALL PRMS (NPRME,IWRK(IPRM+1),MAXP,IHTBL,2*NSTOR) IWRK(IPIN+1) = 0 *------------------------------------------------------------------- * CONTINUE PARTITION THE WORK SPACE *------------------------------------------------------------------- LMAX = (INT(FLOAT(NPRME)*(ALOG10(FLOAT(IWRK(IPRM+1)))- * ALOG10(2.0)) / ALOG10(FLOAT(IPARM(9)))) + 10) * 2 IPIN = IPRM + NPRME IAPRM = IPIN + NPRME KODE = IAPRM + NPRME KNUM = KODE + LMAX KGCD = KNUM + LMAX JWRK = KGCD + LMAX LST = JWRK LENGTH = NWRK - JWRK - 15 LWRK = MAX0(LMAX,LENGTH/5) IW1 = LST + LWRK + 10 IW2 = IW1 + LWRK + 1 IW3 = IW2 + LWRK + 1 IW4 = IW3 + LWRK + 1 *------------------------------------------------------------------- * LOAD THE MATRIX PARAMETERS *------------------------------------------------------------------- REWIND (IPARM(3)) READ (IPARM(3),*) (IWRK(MINFO+I),I=1,4),JBASE IWRK(MINFO+5) = MIN0(IWRK(MINFO+1),IWRK(MINFO+2)) IWRK(MINFO+6) = MAX0(IWRK(MINFO+1),IWRK(MINFO+2)) IWRK(MINFO+7) = 0 IWRK(MINFO+8) = IWRK(MINFO+4) IWRK(MINFO+9) = 0 IWRK(MINFO+10) = IWRK(MINFO+2) + IWRK(MINFO+3) *------------------------------------------------------------------- * INITIALIZE HASH TABLE, ROW AND COLUMN INDEXING, PIVOT VECTOR, * DEPENDENT ROW VECTOR, AND RHS CONSISTENCY VECTOR *------------------------------------------------------------------- DO 60, I = 1, IWRK(MINFO+1) IWRK(IROW+I) = I 60 CONTINUE DO 70, I = 1, IWRK(MINFO+10) IWRK(ICOL+I) = I 70 CONTINUE DO 80, J = 0, NSTOR IHTBL(1,J) = 0 IHTBL(2,J) = 0 80 CONTINUE *------------------------------------------------------------------- * LOAD THE HASH PARAMETERS *------------------------------------------------------------------- IONE = 1 IWRK(IHPRM+1) = NSTOR CALL PRMS (IONE,IP,NSTOR,IWRK(JWRK+1),LENGTH) NSTOR = IP(1) IWRK(IHPRM+2) = 6323 IWRK(IHPRM+3) = IP(1) IWRK(IHPRM+4) = IP(1) IWRK(IHPRM+5) = IWRK(MINFO+2)+IWRK(MINFO+3) IWRK(IHPRM+6) = MAX0 (IWRK(MINFO+1),IWRK(MINFO+2)+1) IWRK(IHPRM+7) = 0 IWRK(IHPRM+8) = 0 IWRK(IHPRM+9) = 0 IWRK(IHPRM+10) = 0 IWRK(IHPRM+11) = ITRM *------------------------------------------------------------------- * LOCAL COMMON VARIABLE INITIALIZATION *------------------------------------------------------------------- ITER = 1 ICPRM = IWRK(IPRM+1) IZERO = 0 IDIAG = 1 IMAT = IPARM(7) IBASE = IPARM(9) IF (IMAT.GT.0) KUN12 = IPARM(4) KUN13 = IPARM(5) ISOLH = 0 ITRI = 0 INULL = 0 IBASIC = 0 JDEN = IWRK(MINFO+2) + 1 IF (IPARM(8).LT.2) ITRI = 1 IF (IPARM(8).EQ.1 .OR. IPARM(8).EQ.3) IBASIC = 1 IF (IPARM(8).EQ.2 .OR. IPARM(8).EQ.3) INULL = 1 II = IFIX(ALOG10(FLOAT(IPARM(9)-1))) + 1 KFORM = (20*IPARM(13)) + II *------------------------------------------------------------------- * PASS # 1 - OPERATIONS * 1. READ DATA MODULO FIRST PRIME * 2. REDUCE TO DIAGONAL FORM WITH PIVOTING * 3. COMPUTE THE HASH TABLE SEPERATION * 4. COMPUTE HADAMARD BOUNDS * 5. STORE RESULTS OF PASS #1 RESULT HASH STRUCTURE *------------------------------------------------------------------- * READ EACH ELEMENT ONE AT A TIME. CONVERT MOD ICPRM AND * STORE IN HASH TABLE. INCREMENT NON-ZERO COUNTER BY COLUMNS. * INSURE THAT ELEMENTS IN THE MATRIX ARE NOT REFERENCED MORE * THAN ONE TIME *------------------------------------------------------------------- IMP = 0 CALL SECOND(TINI) IWRK(IAPRM+1) = ICPRM DO 100, K = 1, IWRK(MINFO+4) READ (IPARM(3),*) I,J,NTRM READ (IPARM(3),*) (IWRK(JWRK+II),II=1,NTRM) IF (I.GT.IWRK(MINFO+1).OR.J.GT.IWRK(IHPRM+5)) THEN WRITE (ITRM,1340) I,J STOP ENDIF IWRK(IPIV+J) = IWRK(IPIV+J) + 1 I1 = MOD(IWRK(JWRK+2),ICPRM) DO 90, II = 3, NTRM I1 = MOD( I1*JBASE + IWRK(JWRK+II),ICPRM) 90 CONTINUE IF (NTRM .GT. 2) IMP = 1 KEY = ISHASH(I,J,0,IWRK(IHPRM+1),IHTBL,NSTOR) IF (IHTBL(2,KEY) .NE. 0) THEN IF ( IER .EQ. 0 ) WRITE (ITRM,1180) WRITE (ITRM,1270) I,J WRITE (ITRM,1190) RETURN ENDIF IHTBL(2,KEY) = I1*IWRK(JWRK+1) 100 CONTINUE *------------------------------------------------------------------- * ROW ELIMINATE TO DIAGONAL FORM USING PIVOTING FROM GJPIV AND * ROW REDUCTION FROM GJRRED. PERFORM INTERMEDIATE OUTPUT IF * SPECIFIED. UPDATE FILLIN INDEX *------------------------------------------------------------------- JZRO = 1 IF (IPARM(7).GE.1) WRITE (ITRM,1060) IWRK(MINFO+1),IWRK(MINFO+2) IF (IPARM(7).GE.2) CALL GJWRT (0,IWRK,NWRK,IHTBL,NSTOR) DO 120, K = 1, IWRK(MINFO+5) KS = K CALL SECOND(T) IF (JZRO .NE. 0) THEN CALL GJPIV (KS,JZRO,IWRK,NWRK,IHTBL,NSTOR) CALL GJRRED (KS,0,IWRK,NWRK,IHTBL,NSTOR) IF (IPARM(7).GE.2) CALL GJWRT (KS,IWRK,NWRK,IHTBL,NSTOR) DO 110, LL = K+1, IWRK(MINFO+2)+IWRK(MINFO+3) KS = KS + IWRK(IPIV+IWRK(ICOL+LL)) 110 CONTINUE IF (KS .GT. IWRK(MINFO+8)) IWRK(MINFO+8) = KS ENDIF 120 CONTINUE IF (IWRK(MINFO+1) .GT. IWRK(MINFO+2)) THEN KK = IWRK(MINFO+6)-IWRK(IDEP+1) DO 125, K = IWRK(MINFO+5)+1, KK, 1 IWRK(IDEP+1) = IWRK(IDEP+1) + 1 IWRK(IDEP+IWRK(IDEP+1)+1) = IWRK(IROW+K) 125 CONTINUE ENDIF IWRK(MINFO+7) = IWRK(MINFO+1) - IWRK(IDEP+1) *------------------------------------------------------------------- * COMPUTE HASH TABLE SPLIT AND RESTORE THE FIRST PASS ELEMENTS IN * ELIMINATION HASH STRUCTURE *------------------------------------------------------------------- KK = INT(FLOAT(IWRK(MINFO+8))/FILL) IF ( IPARM(8) .EQ. 0 .OR. KK .GE. NSTOR ) THEN IPARM(8) = 0 IBASIC = 0 INULL = 0 IPP = IWRK(IHPRM+3) IWRK(IHPRM+12) = NSTOR IWRK(IHPRM+4) = 1 GO TO 160 ENDIF IWRK(IHPRM+14) = 0 DO 130, J = IWRK(MINFO+7)+1, IWRK(MINFO+2)+IWRK(MINFO+3) IWRK(IHPRM+14) = IWRK(IHPRM+14) + IWRK(IPIV+J) 130 CONTINUE RFRAC = FLOAT(IWRK(MINFO+8)) / FLOAT(NSTOR) SFRAC = FLOAT(IWRK(IHPRM+14)*IPARM(12)) / FLOAT(NSTOR*2) TFRAC = RFRAC + SFRAC IF ( RFRAC .GT. .95 ) THEN MAXP = MIN0(INT(FLOAT(IWRK(MINFO+8))/.9),IWRK(IHPRM+3)) ELSEIF ( TFRAC .GE. .9 ) THEN MAXP = MIN0(INT(FLOAT(IWRK(MINFO+8))/FILL),IWRK(IHPRM+3)) ELSE TFRAC = (1.-TFRAC)*(2./ 3.) SFRAC = FLOAT(NSTOR) MAXP = INT((TFRAC+RFRAC)*SFRAC) ENDIF CALL PRMS (IONE,IP,MAXP,IWRK(JWRK+1),LENGTH) IWRK(IHPRM+12) = IP(1) IPP = IWRK(IHPRM+3) DO 150, I = 1, IWRK(MINFO+1) DO 140, J = I, IWRK(MINFO+2)+IWRK(MINFO+3) IWRK(IHPRM+3) = IPP KEY = IRHASH(IWRK(IROW+I),IWRK(ICOL+J),0,IWRK(IHPRM+1), * IHTBL,NSTOR) KK = IHTBL(2,KEY) IF ( KK .NE. 0 ) THEN IHTBL(2,KEY) = 0 IWRK(IHPRM+3) = IP(1) KEY = ISHASH(IWRK(IROW+I),IWRK(ICOL+J),0,IWRK(IHPRM+1), * IHTBL,NSTOR) IHTBL(2,KEY) = KK IF ( I.GT.IWRK(MINFO+7) .AND. J.GT.IWRK(MINFO+2) ) THEN I1 = J - IWRK(MINFO+2) IWRK(IRHS+I1) = 1 ENDIF ENDIF 140 CONTINUE 150 CONTINUE IWRK(IHPRM+3) = IPP IPP = IP(1) MHZ = IP(1) + 1 MAXP = MAX(NSTOR-IP(1)-1,1) IP(1) = 0 CALL PRMS (IONE,IP,MAXP,IWRK(JWRK+1),LENGTH) IWRK(IHPRM+4) = IP(1) 160 CONTINUE *------------------------------------------------------------------- * COMPUTE HADAMARD BOUNDS IF NECESSARY. READ EACH ELEMENT ONE AT * A TIME. COMPUTE REGULAR HADAMARD BOUND. KEEP ROW AND COLUMN * NORMS, AND KEEP THE RIGHT HAND SIDE NORMS. *------------------------------------------------------------------- IF ( IPARM(10).LE.0 ) THEN DO 165, I = 1, NH HWRK(I) = 0. 165 CONTINUE REWIND (IPARM(3)) READ (IPARM(3),*) (IWRK(JWRK+II),II=1,4),JBASE HBASE = FLOAT(JBASE) HLOG = ALOG10(HBASE) DO 170, K = 1, IWRK(MINFO+4) READ (IPARM(3),*) I,J,NTRM READ (IPARM(3),*) (IWRK(JWRK+II),II=1,NTRM) NTRM = (NTRM-2)*2 ICOEF = (IWRK(JWRK+2)+IMP)**2 IF ( IPARM(10) .NE. 0 ) THEN KEY = ISHASH(I,J,1,IWRK(IHPRM+1),IHTBL,NSTOR) IHTBL(2,KEY) = NTRM KEY = ISHASH(I,J,2,IWRK(IHPRM+1),IHTBL,NSTOR) IHTBL(2,KEY) = ICOEF ENDIF *------------------------------------------------------------------- * INCREMENT NORMS WITH PSEUDO-FLOATING POINT ARITHMETIC BY STORING * THE COEFFICIENT AND THE POWER OF THE ROW/COLUMN NORMS IN A * REAL AND INTEGER ARRAY RESPECTIVELY. *------------------------------------------------------------------- IF (J .GT. IWRK(MINFO+2)) THEN J = J + IW3 + IWRK(MINFO+1) JJ = J - IW3 IDIF = IABS(IWRK(J) - NTRM) IF (IWRK(J).LT.NTRM) THEN IDIF = MAX0(-3,IWRK(J)-NTRM) HWRK(JJ) = FLOAT(ICOEF) + (HWRK(JJ)*(HBASE**IDIF)) IWRK(J) = NTRM ELSE IDIF = MAX0(-3,NTRM-IWRK(J)) HWRK(JJ) = HWRK(JJ) + (FLOAT(ICOEF)*(HBASE**IDIF)) ENDIF ELSE J = J + IWRK(MINFO+1) IDIF = IABS(IWRK(IW3+J) - NTRM) IF ( IWRK(IW3+J).LT.NTRM ) THEN IDIF = MAX0(-3,IWRK(IW3+J)-NTRM) HWRK(J) = FLOAT(ICOEF) + (HWRK(J)*(HBASE**IDIF)) IWRK(IW3+J) = NTRM ELSE IDIF = MAX0(-3,NTRM-IWRK(IW3+J)) HWRK(J) = HWRK(J) + (FLOAT(ICOEF)*(HBASE**IDIF)) ENDIF IF (IWRK(IW3+I) .LT. NTRM ) THEN IDIF = MAX0(-3,IWRK(IW3+I)-NTRM) HWRK(I) = FLOAT(ICOEF) + (HWRK(I)*(HBASE**IDIF)) IWRK(IW3+I) = NTRM ELSE IDIF = MAX0(-3,NTRM-IWRK(IW3+I)) HWRK(I) = HWRK(I) + (FLOAT(ICOEF)*(HBASE**IDIF)) ENDIF ENDIF 170 CONTINUE *------------------------------------------------------------------- * COMPUTE BOUND OF HADAMARD. COMPUTE THE MAX. COLUMN NORM OF * THE RIGHT HAND SIDES. *------------------------------------------------------------------- HDMRD = 1.0 IMAX = 0 DO 180, J = IWRK(MINFO+2) + 1, IWRK(MINFO+2)+IWRK(MINFO+3) JJ = J + IW3 + IWRK(MINFO+1) KK = JJ - IW3 IF (IWRK(JJ) .GT. IMAX) THEN IMAX = IWRK(JJ) HDMRD = HWRK(KK) ELSEIF (IWRK(JJ) .EQ. IMAX) THEN HDMRD = MAX(HWRK(KK),HDMRD) ENDIF 180 CONTINUE HDMRD = MAX(HDMRD,1.) BMAX = (ALOG10(2.*HDMRD) + (FLOAT(IMAX)*HLOG)) / 2. *------------------------------------------------------------------- * BOUND THE DETERMINANT OF THE ROW NORMS IN THE MATRIX. *------------------------------------------------------------------- HDMRD = BMAX DO 190, I = 1, IWRK(MINFO+1) IF ( IWRK(MINFO+2) .GE. IWRK(MINFO+1) ) THEN IMAX = IWRK(IW3+I) HMAX = MAX(1.,HWRK(I)) LL = I ELSE IMAX = 0 HMAX = 1. DO 185, J = 1, IWRK(MINFO+1) IF ( IWRK(IW3+J) .GT. IMAX ) THEN IMAX = IWRK(IW3+J) HMAX = HWRK(J) LL = J ELSEIF ( (IWRK(IW3+J).EQ.IMAX) .AND. * (HWRK(J).GT.HMAX)) THEN HMAX = HWRK(J) LL = J ENDIF 185 CONTINUE ENDIF HMAX = MAX(HMAX,1.) HDMRD = HDMRD + ((ALOG10(HMAX) + (FLOAT(IMAX)*HLOG)) / 2.) IWRK(IW3+LL) = 0 HWRK(LL) = 0. 190 CONTINUE *------------------------------------------------------------------- * COMPUTE NUMBER OF PRIMES REQUIRED TO BOUND HDMRD BY ROWS *------------------------------------------------------------------- HMAX = 0. DO 195, I = 1, NPRME PRIME = ALOG10(FLOAT(IWRK(IPRM+I))) HMAX = HMAX + PRIME IF ( HDMRD.LE.HMAX ) THEN NPASS = I GO TO 200 ENDIF 195 CONTINUE NPASS = NPRME 200 CONTINUE *------------------------------------------------------------------- * BOUND THE DETERMINANT OF THE COLUMN NORMS IN THE MATRIX. *------------------------------------------------------------------- HDMRD = BMAX DO 210, J = IWRK(MINFO+1)+1, IWRK(MINFO+1)+IWRK(MINFO+2) IF ( IWRK(MINFO+1) .GE. IWRK(MINFO+2) ) THEN IMAX = IWRK(IW3+J) HMAX = MAX(1.,HWRK(J)) ELSE IMAX = 0 HMAX = 1. DO 205, I=IWRK(MINFO+1)+1, IWRK(MINFO+1)+IWRK(MINFO+2) IF ( IWRK(IW3+I) .GT. IMAX ) THEN IMAX = IWRK(IW3+I) HMAX = HWRK(I) LL = I ELSEIF ( (IWRK(IW3+I).EQ.IMAX) .AND. * (HWRK(I).GT.HMAX)) THEN HMAX = HWRK(I) LL = I ENDIF 205 CONTINUE ENDIF HMAX = MAX(HMAX,1.) HDMRD = HDMRD + ((ALOG10(HMAX) + (FLOAT(IMAX)*HLOG)) / 2.) IWRK(IW3+LL) = 0 HWRK(LL) = 0. 210 CONTINUE *------------------------------------------------------------------- * COMPUTE NUMBER OF PRIMES REQUIRED TO BOUND HDMRD BY COLUMNS *------------------------------------------------------------------- HMAX = 0. DO 215, I = 1, NPRME PRIME = ALOG10(FLOAT(IWRK(IPRM+I))) HMAX = HMAX + PRIME IF ( HDMRD.LE.HMAX ) THEN NPASS1 = I GO TO 220 ENDIF 215 CONTINUE NPASS1 = NPRME 220 CONTINUE *------------------------------------------------------------------- * TAKE THE SMALLEST NUMBER OF PRIMES THAT BOUND THE NORMS *------------------------------------------------------------------- NPASS = MIN0(NPASS,NPASS1) NPASS1 = 0 IF ( NPASS .GE. NPRME ) NPASS2 = 1 ENDIF *------------------------------------------------------------------- * COMPUTE THE MODIFIED HADAMARD BOUND ON THE RXR SQUARE SUBMATRIX * THAT IS NONSINGULAR AND CHOSEN BY PIVOTING. WORK THROUGH THE * INDEPENDENT ROWS AND COLUMNS TO COMPUTE ROW NORMS OF THE * NON-SINGULAR SUB-MATRIX. INCREMENT HADAMARD BOUND, HDMRD. *------------------------------------------------------------------- IF ( IPARM(10).LE.-1 .AND. IWRK(MINFO+7).EQ.IWRK(MINFO+1) ) THEN NPASS1 = NPASS ELSEIF ( IPARM(10) .LE. -1 ) THEN HDMRD = BMAX DO 240, I = 1, IWRK(MINFO+7) HMAX = 1. IMAX = 0 DO 230, J = 1, IWRK(MINFO+7) KEY = IRHASH(IWRK(IROW+I),IWRK(ICOL+J),1,IWRK(IHPRM+1), * IHTBL,NSTOR) NTRM = IHTBL(2,KEY) KEY = IRHASH(IWRK(IROW+I),IWRK(ICOL+J),2,IWRK(IHPRM+1), * IHTBL,NSTOR) HINT = FLOAT(IHTBL(2,KEY)) IF ( NTRM .EQ. IMAX ) THEN HMAX = HMAX + HINT ELSEIF ( NTRM .GT. IMAX ) THEN IDIF = IMAX - NTRM HMAX = HINT + (HMAX*(HBASE**IDIF)) IMAX = NTRM ELSE IDIF = NTRM - IMAX HMAX = HMAX + (HINT*(HBASE**IDIF)) ENDIF 230 CONTINUE HMAX = MAX(HMAX,1.) HDMRD = HDMRD + ((ALOG10(HMAX)+(FLOAT(IMAX)*HLOG)) / 2.) 240 CONTINUE *------------------------------------------------------------------- * COMPUTE NUMBER OF PRIMES REQUIRED TO BOUND MODIFIED HDMRD *------------------------------------------------------------------- HMAX = 0. DO 250, I = 1, NPRME PRIME = ALOG10(FLOAT(IWRK(IPRM+I))) HMAX = HMAX + PRIME IF ( HDMRD.LE.HMAX ) THEN NPASS1 = I GO TO 260 ENDIF 250 CONTINUE NPASS1 = NPRME NPASS3 = 1 260 CONTINUE ENDIF *------------------------------------------------------------------- * COMPUTE INITIAL RESULTS AND ZERO THE HADAMARD PART OF HASHING *------------------------------------------------------------------- IWRK(IHPRM+3) = IPP II = IWRK(IHPRM+5)*IWRK(IHPRM+6) DO 270, I= 0, NSTOR IF (IHTBL(1,I).GT.II .OR. I.GT.IPP) THEN IHTBL(1,I) = 0 IHTBL(2,I) = 0 ENDIF 270 CONTINUE CALL GJRES (IWRK,NWRK,IHTBL,NSTOR) *------------------------------------------------------------------- * COMPUTE AND OUTPUT INITIAL RESULTS *------------------------------------------------------------------- IF (IWRK(MINFO+7).EQ.IWRK(MINFO+5) .AND. IPARM(7).GE.1) THEN WRITE (ITRM,1030) IF ( IPARM(10) .LE. 0 ) WRITE (ITRM,1010) NPASS IF ( NPASS2 .EQ. 1 ) WRITE (ITRM,1025) IF (IWRK(MINFO+1).GT.IWRK(MINFO+2)) THEN WRITE (ITRM,1045) IWRK(IDEP+1) WRITE (ITRM,9080) (IWRK(IDEP+I),I=2,IWRK(IDEP+1)+1) ENDIF WRITE (ITRM,9090) ELSEIF (IPARM(7).GE.1 .AND. IPARM(10).LE.0) THEN WRITE (ITRM,1000) WRITE (ITRM,1010) NPASS IF ( NPASS1 .GT. 0 ) WRITE (ITRM,1020) NPASS1 IF ( NPASS2 .EQ. 1 ) WRITE (ITRM,1025) IF ( NPASS3 .EQ. 1 ) WRITE (ITRM,1035) WRITE (ITRM,1040) IWRK(IDEP+1) WRITE (ITRM,9080) (IWRK(IDEP+I),I=2,IWRK(IDEP+1)+1) WRITE (ITRM,9090) ELSEIF (IPARM(7).GE.1) THEN WRITE (ITRM,1000) WRITE (ITRM,1040) IWRK(IDEP+1) WRITE (ITRM,9080) (IWRK(IDEP+I),I=2,IWRK(IDEP+1)+1) WRITE (ITRM,9090) ENDIF *------------------------------------------------------------------- * CALCULATE THE REQUIRED NUMBER OF PASSES *------------------------------------------------------------------- IF ((IPARM(8).EQ.0) .AND. (IWRK(MINFO+7).EQ.IWRK(MINFO+1)) .AND. * (IWRK(MINFO+2).GE.IWRK(MINFO+1))) THEN ITER = 2 KOUNT = 0 GO TO 360 ELSEIF (IPARM(10) .GT. 1) THEN NTERM = IPARM(10) ELSEIF (IPARM(10) .EQ. 1) THEN ITER = 2 KOUNT = 0 GO TO 360 ELSEIF (IPARM(10) .EQ. 0) THEN NTERM = NPASS ELSEIF (IPARM(10) .EQ. -1) THEN NTERM = NPASS ELSEIF (IPARM(10) .EQ. -2) THEN NTERM = NPASS1 ENDIF *------------------------------------------------------------------- * PASS # 2,3, ... - OPERATIONS BEGIN WITH * 1. SET PRIMES * 2. INITIALIZE DATA STRUCTURES * 3. READ DATA MODULO ITERTH PRIME * 4. REDUCE TO IDENTITY FORM -- NO PIVOTING * 5. WRITE PASS ITER RESULTS TO INTERMEDIATE FILES *------------------------------------------------------------------- NDIS = 0 KOUNT = 0 DO 350, ITER = 2, NTERM 280 CONTINUE ICPRM = IWRK(IPRM+ITER+NDIS) IWRK(IAPRM+ITER) = ICPRM IPP = 1 DO 290, K = 1, ITER-1 IPP = MOD(IPP*IWRK(IAPRM+K),ICPRM) 290 CONTINUE IWRK(IPIN+ITER) = INV(IPP,ICPRM) DO 300, I = 0, IWRK(IHPRM+3) IHTBL(1,I) = 0 IHTBL(2,I) = 0 300 CONTINUE DO 310, I = 1, IWRK(MINFO+1) IWRK(IPIV+I) = 0 310 CONTINUE *------------------------------------------------------------------- * READ EACH ELEMENT ONE AT A TIME. CONVERT MOD ICPRM AND * STORE IN HASH TABLE. INCREMENT NON-ZERO COUNTER BY ROWS. *------------------------------------------------------------------- REWIND (IPARM(3)) READ (IPARM(3),*) (IWRK(JWRK+II),II=1,4),JBASE DO 330, K = 1, IWRK(MINFO+4) READ (IPARM(3),*) I,J,NTRM READ (IPARM(3),*) (IWRK(JWRK+II),II=1,NTRM) IWRK(IPIV+I) = IWRK(IPIV+I) + 1 I1 = MOD(IWRK(JWRK+2),ICPRM) DO 320, L = 3, NTRM I1 = MOD( I1*JBASE + IWRK(JWRK+L),ICPRM) 320 CONTINUE KEY = ISHASH(I,J,0,IWRK(IHPRM+1),IHTBL,NSTOR) IHTBL(2,KEY) = I1*IWRK(JWRK+1) 330 CONTINUE *------------------------------------------------------------------- * ELIMINATE TO APPROPRIATE FORM - PERFORM OUTPUT IF DESIRED *------------------------------------------------------------------- IDIAG = 1 IF (IPARM(7).GE.4) CALL GJWRT (0,IWRK,NWRK,IHTBL,NSTOR) DO 340, K = 1, IWRK(MINFO+7) KS = K CALL GJRRED (KS,1,IWRK,NWRK,IHTBL,NSTOR) IF (IPARM(7).GE.4) CALL GJWRT (KS,IWRK,NWRK,IHTBL,NSTOR) 340 CONTINUE *------------------------------------------------------------------- * STORE NEW RESIDUES AND PERFORM RANK CHECK - * IF IRANK = 1 RANK INCREASED - DISCARD ALL PRIMES * IF IRANK = -1 RANK DECREASED - DISCARD THIS PRIME * IF IRANK = 0 RANK UNCHANGED - CONTINUE *------------------------------------------------------------------- CALL GJRES (IWRK,NWRK,IHTBL,NSTOR) IF ( ISOLH .EQ. 1 ) THEN IF ( IPARM(7) .GE. 1 ) WRITE (ITRM,1095) IBASIC = 0 INULL = 0 KOUNT = 0 GO TO 360 ENDIF IF ( IRANK .EQ. 1 ) THEN IF ( IPARM(7) .GE. 1 ) WRITE (ITRM,1080) ICPRM NDIS = 0 MAXP = ICPRM-2 GO TO 30 ELSEIF ( IRANK .EQ. -1 ) THEN IF ( IPARM(7) .GE. 1 ) WRITE (ITRM,1090) ICPRM NDIS = NDIS + 1 GO TO 280 ENDIF IF ( IPARM(7) .GE. 1 ) THEN CALL SECOND(T) WRITE (ITRM,1050) ITER,T-TINI ENDIF *------------------------------------------------------------------- * PERFORM CONVERGENCE CHECK: STOP IF * 1. THREE SUCCESSIVE ITERATIONS DO NOT CHANGE THE RESULTS * 2. THE MODIFIED HADAMARD BOUND IS SATISFIED AND THE RESULTS * DONT CHANGE ON ONE ITERATE. *------------------------------------------------------------------- IF ( IZERO .NE. 0 ) THEN KOUNT = 0 ELSE KOUNT = KOUNT + 1 ENDIF IF ( IPARM(10).EQ.-1 .AND. (KOUNT.EQ.4 .OR. * (KOUNT.EQ.1 .AND. ITER.GT.NPASS1)) ) GO TO 360 350 CONTINUE *------------------------------------------------------------------- * BEGIN OUTPUT - WRITE HEADER AND MATRIX DEPENDENCIES *------------------------------------------------------------------- 360 CONTINUE ITTOL = ITER - KOUNT IF ( KUN12 .GT. 0 ) THEN CALL GOLIA(KUN12) WRITE (KUN12,9000) WRITE (KUN12,9010) IWRK(MINFO+1),IWRK(MINFO+2),IWRK(MINFO+3) IF (NPASS.GT.0) WRITE (KUN12,9015) NPASS IF (NPASS1.GT.0) WRITE (KUN12,9020) NPASS1 WRITE (KUN12,9025) ITER-1 WRITE (KUN12,9045) KOUNT IF ( IWRK(MINFO+7) .EQ. IWRK(MINFO+1) ) THEN WRITE (KUN12,9035) IWRK(MINFO+7) ELSE WRITE (KUN12,9030) IWRK(MINFO+7) WRITE (KUN12,9080) (IWRK(IDEP+I),I=2,IWRK(IDEP+1)+1) ENDIF WRITE (KUN12,9040) WRITE (KUN12,9080) (IWRK(IROW+I),I=1,IWRK(MINFO+7)) WRITE (KUN12,9050) WRITE (KUN12,9080) (IWRK(ICOL+I),I=1,IWRK(MINFO+7)) ENDIF *------------------------------------------------------------------- * COMPUTE AND OUTPUT THE MINOR OF THE LEADING SQUARE SYSTEM *------------------------------------------------------------------- IF (KUN12.GT.0 .AND. IPARM(13).EQ.1) WRITE (KUN12,9060) IBASE IF (KUN12.GT.0 .AND. IPARM(13).EQ.0) WRITE (KUN12,9065) IF ( KUN12.GT.0 .AND. IWRK(MINFO+7).EQ.IWRK(MINFO+1) .AND. * IWRK(MINFO+2).EQ.IWRK(MINFO+1) ) THEN WRITE (KUN12,9075) ELSEIF ( KUN12 .GT. 0 ) THEN WRITE (KUN12,9070) ENDIF CALL RADIXF (KNUM,ITTOL,KODE,IWRK,NWRK) IF ( TOUT ) CALL MPOUT (IWRK(KODE+1),LMAX,KUN12,KFORM) *------------------------------------------------------------------- * COMPUTE AND OUTPUT THE DETERMINANT OF THE RATIONAL MATRIX * IF THE INTEGER MATRIX FILE CONTAINS THE DETERMINANT OF THE * DIAGONAL SCALING MATRIX USED CONVERT IT. *------------------------------------------------------------------- READ (IPARM(3),*,END=370) I,J,NTRM IF ( (I.EQ.0) .AND. (J.EQ.0) .AND. (KUN12.GT.0) .AND. * (LWRK.GT.NTRM)) THEN WRITE (KUN12,9085) READ (IPARM(3),*) (IWRK(JWRK+II),II=0,NTRM-1) IF (IBASE.NE.JBASE) * CALL MPBASE (IWRK(JWRK),NTRM,JBASE,IBASE,IWRK(IW2)) IWRK(JWRK) = IWRK(JWRK)*NTRM IF (JBASE.NE.IBASE) * CALL MPBASE (IWRK(JWRK),NTRM+1,IBASE,JBASE,IWRK(IW3)) IF ( LWRK.LE.NTRM .OR. IPARM(10).EQ.1 ) THEN CALL MPOUT (IWRK(KODE+1),LMAX,KUN12,KFORM) WRITE (KUN12,9270) CALL MPOUT (IWRK(JWRK),LENGTH,KUN12,KFORM) ELSE CALL MPCOPY (IWRK(KGCD+1),IWRK(JWRK),LMAX) CALL MPCOPY (IWRK(IW4),IWRK(KGCD+1),LMAX) CALL MPCOPY (IWRK(KNUM+1),IWRK(KODE+1),LMAX) CALL MPGCD (IWRK(KNUM+1),IWRK(KGCD+1),LMAX,IBASE, * IWRK(IW1),IWRK(IW2),IWRK(IW3),IWRK(LST)) CALL MPDIV (IWRK(KODE+1),IWRK(KGCD+1),IWRK(LST+1), * IWRK(IW1),LMAX,IBASE,IWRK(IW2),IWRK(IW3)) CALL MPOUT (IWRK(LST+1),LMAX,KUN12,KFORM) CALL MPDIV (IWRK(IW4),IWRK(KGCD+1),IWRK(LST+1), * IWRK(IW1),LMAX,IBASE,IWRK(IW2),IWRK(IW3)) WRITE (KUN12,9270) CALL MPOUT (IWRK(LST+1),LMAX,KUN12,KFORM) ENDIF ENDIF 370 CONTINUE *------------------------------------------------------------------- * ZERO ELIMINATION HASH TABLE AND RESTRUCTURE FOR OUTPUT *------------------------------------------------------------------- IWRK(IHPRM+13) = 0 DO 380, I = 0, IWRK(IHPRM+3) IHTBL(1,I) = 0 IHTBL(2,I) = 0 380 CONTINUE DO 390, I = MHZ, MHZ+IWRK(IHPRM+4) IF (IHTBL(1,I).GT.0) IWRK(IHPRM+13) = IWRK(IHPRM+13) + 1 390 CONTINUE IWRK(IHPRM+3) = NSTOR *------------------------------------------------------------------- * COMPUTE AND OUTPUT PARTICULAR SOLUTIONS. TEST EACH RHS FOR * CONSISTENTCY. *------------------------------------------------------------------- IWRK(IHPRM+11) = 0 IF (IBASIC .EQ. 1 ) THEN IF ( TOUT ) WRITE (KUN12,9200) IF ( IPARM(7) .GE. 1 ) WRITE (ITRM,1300) DO 470, J = 1, IWRK(MINFO+3) JSOL = IWRK(MINFO+2) + J IF ( IWRK(IRHS+J) .EQ. 1 ) THEN IF ( TOUT ) WRITE (KUN12,9210) J IF ( IPARM(7) .GE. 1 ) WRITE (ITRM,1310) J GO TO 470 ENDIF *------------------------------------------------------------------- * IF JTH SOLUTION VECTOR ELEMENTS ARE NON-ZERO, FIND THE RESIDUE * MODULO ICPRM AND COMPUTE THE NEW MIXED RADIX REPRESENTATION. * UPDATE THE GCD OF THE SOLUTIONS NUMERATORS AND DENOMINATOR. * RESTORE SOLUTION ELEMENTS WHILE GCD IS BEING COMPUTED. *------------------------------------------------------------------- IF ( TOUT ) WRITE (KUN12,9220) J, JSOL CALL MPCOPY(IWRK(KGCD+1),IWRK(KODE+1),LMAX) DO 410, I = 1, IWRK(MINFO+2) IS = I CALL JRHASH (IS,JSOL,LST,ITTOL,0,JZ,IWRK,NWRK, * IHTBL,NSTOR) IF ( JZ .EQ. 0 ) THEN IWRK(IPIV+I) = 0 ELSE IWRK(IPIV+I) = 1 CALL RADIXF (LST,ITTOL,KNUM,IWRK,NWRK) DO 400, K = 1, IABS(IWRK(KNUM+1)), 1 KS = K KEY = ISHASH(IS,JSOL,KS,IWRK(IHPRM+1),IHTBL,NSTOR) IHTBL(2,KEY) = IWRK(KNUM+K) 400 CONTINUE CALL MPGCD (IWRK(KNUM+1),IWRK(KGCD+1),LMAX,IBASE, * IWRK(IW1),IWRK(IW2),IWRK(IW3),IWRK(LST)) ENDIF 410 CONTINUE *------------------------------------------------------------------- * COMPUTE AND OUTPUT THE GCD FOR JTH RHS. STORE DENOMINATOR * FOR CURRENT RIGHT HAND SIDE IN HASH TABLE. *------------------------------------------------------------------- IF ( TOUT ) WRITE (KUN12,9270) CALL MPDIV (IWRK(KODE+1),IWRK(KGCD+1),IWRK(IW1), * IWRK(IW2),LMAX,IBASE,IWRK(IW3),IWRK(IW4)) IF ( TOUT ) CALL MPOUT (IWRK(IW1),LMAX,KUN12,KFORM) DO 420, K = 0, IABS(IWRK(IW1))-1, 1 KS = K+1 KEY = ISHASH(JDEN,JSOL,KS,IWRK(IHPRM+1),IHTBL,NSTOR) IHTBL(2,KEY) = IWRK(IW1+K) 420 CONTINUE *------------------------------------------------------------------- * COMPUTE AND OUTPUT THE NUMERATORS BY DIVIDING OUT THE GCD. * STORE NUMERATORS IN THE HASH TABLE. *------------------------------------------------------------------- IF ( TOUT ) WRITE (KUN12,9260) DO 460, I = 1, IWRK(MINFO+2) IS = I IF ( IWRK(IPIV+IS) .NE. 0 ) THEN KEY = IRHASH(IS,JSOL,1,IWRK(IHPRM+1),IHTBL,NSTOR) IWRK(LST+1) = IHTBL(2,KEY) IHTBL(2,KEY) = 0 DO 440, K = 2, IABS(IWRK(LST+1)), 1 KS = K KEY = IRHASH(IS,JSOL,KS,IWRK(IHPRM+1),IHTBL,NSTOR) IWRK(LST+K) = IHTBL(2,KEY) IHTBL(2,KEY) = 0 440 CONTINUE CALL MPDIV (IWRK(LST+1),IWRK(KGCD+1),IWRK(IW1), * IWRK(IW2),LMAX,IBASE,IWRK(IW3),IWRK(IW4)) IF ( TOUT ) THEN WRITE (KUN12,9140) IS CALL MPOUT (IWRK(IW1),LMAX,KUN12,KFORM) END IF DO 450, K = 0, IABS(IWRK(IW1))-1, 1 KS = K+1 KEY = ISHASH(IS,JSOL,KS,IWRK(IHPRM+1),IHTBL,NSTOR) IHTBL(2,KEY) = IWRK(IW1+K) 450 CONTINUE ENDIF 460 CONTINUE *------------------------------------------------------------------- * CALCULATE TIME AND GO TO NEXT NUMERATOR. *------------------------------------------------------------------- IF ( IPARM(7) .GE. 1 ) THEN CALL SECOND(T) WRITE (ITRM,1310) J,T-TINI ENDIF 470 CONTINUE ENDIF *------------------------------------------------------------------- * COMPUTE AND OUTPUT NULL SPACE *------------------------------------------------------------------- IF ( IWRK(MINFO+7).EQ.IWRK(MINFO+2) .AND. INULL.EQ.1 ) THEN IF ( TOUT ) WRITE (KUN12,9110) ELSEIF ( IWRK(MINFO+7).LT.IWRK(MINFO+2) .AND. INULL.EQ.1 ) THEN IF ( IPARM(7) .GE. 1 ) WRITE (ITRM,1320) * IWRK(MINFO+2)-IWRK(MINFO+7) IF ( TOUT ) WRITE (KUN12,9100) KK = IWRK(MINFO+2) - IWRK(MINFO+7) IF ( TOUT ) WRITE (KUN12,9130) KK IF ( TOUT ) WRITE (KUN12,9080) * (IWRK(ICOL+I), I = IWRK(MINFO+7)+1, IWRK(MINFO+2)) *------------------------------------------------------------------- * REDUCE ELEMENTS IN EACH VECTOR IN NULL SPACE TO LOWEST TERMS *------------------------------------------------------------------- DO 540, J = 1, IWRK(MINFO+2)-IWRK(MINFO+7) JS = J JSOL = J + IWRK(MINFO+7) LSOL = IWRK(ICOL+JSOL) CALL MPCOPY (IWRK(KGCD+1),IWRK(KODE+1),LMAX) IF ( TOUT ) WRITE (KUN12,9120) J,JSOL,IWRK(ICOL+JSOL) *------------------------------------------------------------------- * IN JTH VECTOR, ELEMENT-WISE COMPUTE THE FIXED RADIX FORM. * UPDATE THE GCD OF NULL VECTOR ELEMENTS. RESTORE EACH ELEMENT. *------------------------------------------------------------------- DO 500, I = 1, IWRK(MINFO+2) IS = I IF ( I .EQ. LSOL ) THEN CALL MPCOPY (IWRK(KNUM+1),IWRK(KODE+1),LMAX) IWRK(KNUM+1) = (-1)*IWRK(KNUM+1) ELSE CALL JRHASH (IS,JSOL,LST,ITTOL,0,IZ,IWRK,NWRK, * IHTBL,NSTOR) IF ( IZ .EQ. 0 ) GO TO 500 CALL RADIXF (LST,ITTOL,KNUM,IWRK,NWRK) ENDIF DO 480, K = 1, IABS(IWRK(KNUM+1)), 1 KS = K KEY = ISHASH(IS,JSOL,KS,IWRK(IHPRM+1),IHTBL,NSTOR) IHTBL(2,KEY) = IWRK(KNUM+K) 480 CONTINUE CALL MPGCD (IWRK(KNUM+1),IWRK(KGCD+1),LMAX,IBASE, * IWRK(IW1),IWRK(IW2),IWRK(IW3),IWRK(LST)) 500 CONTINUE *------------------------------------------------------------------- * OUTPUT ELEMENTS OF THE JTH VECTOR AFTER DIVIDING OUT THE GCD. * STORE IN HASH TABLE. *------------------------------------------------------------------- DO 530, I = 1, IWRK(MINFO+2) IS = I KEY = IRHASH(IS,JSOL,1,IWRK(IHPRM+1),IHTBL,NSTOR) IWRK(LST+1) = IHTBL(2,KEY) IHTBL(2,KEY) = 0 DO 510, K = 2, IABS(IWRK(LST+1)), 1 KS = K KEY = IRHASH(IS,JSOL,KS,IWRK(IHPRM+1),IHTBL,NSTOR) IWRK(LST+K) = IHTBL(2,KEY) IHTBL(2,KEY) = 0 510 CONTINUE IF ( IWRK(LST+1) .NE. 0 ) THEN CALL MPDIV (IWRK(LST+1),IWRK(KGCD+1),IWRK(IW1), * IWRK(IW2),LMAX,IBASE,IWRK(IW3),IWRK(IW4)) IF ( TOUT ) WRITE (KUN12,9140) IS IF ( TOUT ) CALL MPOUT (IWRK(IW1),LMAX,KUN12,KFORM) DO 520, K = 0, IABS(IWRK(IW1))-1, 1 KS = K+1 KEY = ISHASH(IS,JSOL,KS,IWRK(IHPRM+1),IHTBL,NSTOR) IHTBL(2,KEY) = IWRK(IW1+K) 520 CONTINUE ENDIF 530 CONTINUE *------------------------------------------------------------------- * CALCULATE TIME AND GO TO NEXT ELEMENT. *------------------------------------------------------------------- IF ( IPARM(7) .GE. 1 ) THEN CALL SECOND(T) WRITE (ITRM,1330) J,T-TINI ENDIF 540 CONTINUE ENDIF *------------------------------------------------------------------- * COMPUTE AND OUTPUT HASHING STATS *------------------------------------------------------------------- IF ( KUN12 .GT. 0 ) THEN WRITE (KUN12,9300) WRITE (KUN12,9330) IWRK(IHPRM+7) HWRK(1) = (FLOAT(IWRK(MINFO+8))*100.) / FLOAT(IWRK(IHPRM+12)) WRITE (KUN12,9340) IWRK(MINFO+8),IWRK(IHPRM+12),HWRK(1) ITBL = IWRK(MINFO+8)+100 IF ( IPARM(8) .NE. 0 ) THEN HWRK(1) = FLOAT(IWRK(IHPRM+13)*100) / FLOAT(IWRK(IHPRM+4)) WRITE (KUN12,9345) IWRK(IHPRM+13),IWRK(IHPRM+4),HWRK(1) WRITE (KUN12,9335) IWRK(IHPRM+8) ITBL = ITBL+IWRK(IHPRM+13) ENDIF IF (IWRK(IHPRM+11) .EQ. -1) WRITE (KUN12,9355) HWRK(1) = FLOAT(IWRK(IHPRM+10)+IWRK(IHPRM+9)) / * FLOAT(IWRK(IHPRM+10)) WRITE (KUN12,9310) HWRK(1) HWRK(1) = FLOAT(IWRK(MINFO+8))/FLOAT(IWRK(MINFO+4)) WRITE (KUN12,9320) HWRK(1) HWRK(1) = FLOAT(IWRK(MINFO+9)+IWRK(MINFO+4)) / * FLOAT(IWRK(MINFO+4)) WRITE (KUN12,9350) HWRK(1) ENDIF *------------------------------------------------------------------- * CLOSE DOWN - STORE THE MINOR IN THE HASH TABLE, LOAD IPARM * ARRAY WITH HASHING PARAMETERS AND RETURN. *------------------------------------------------------------------- DO 560, K = 1, IABS(IWRK(KODE+1)), 1 KS = K KEY = ISHASH(JDEN,1,KS,IWRK(IHPRM+1),IHTBL,NSTOR) IHTBL(2,KEY) = IWRK(KODE+K) 560 CONTINUE IF (IWRK(IHPRM+11) .EQ. -1) WRITE (ITRM,9355) IPARM(0) = 0 IF ( ISOLH .EQ. 1 ) IPARM(0) = -2 DO 570, K = 1, 11 IPARM(K) = IWRK(IHPRM+K) 570 CONTINUE IPARM(4) = IWRK(MINFO+7) KK = 0 DO 580, K = 0, IWRK(IHPRM+3)-1 IF ( IHTBL(2,K) .NE. 0 ) KK = KK + 1 580 CONTINUE HWRK(1) = FLOAT(KK)*100. / FLOAT(IWRK(IHPRM+3)) IF ( TOUT ) WRITE (KUN12,9360) KK,IWRK(IHPRM+3),HWRK(1) IF ( KUN12 .GT. 0 .AND. ISOLH .EQ. 0 ) THEN ITBL = INT(FLOAT(ITBL+KK)/FILL)*2 + NWRK WRITE (KUN12,9370) ITBL ELSEIF ( KUN12 .GT. 0 ) THEN ITBL = (2*IWRK(MINFO+8)) + (IWRK(MINFO+14)*NTERM) + NWRK WRITE (KUN12,9380) ITBL ENDIF RETURN *------------------------------------------------------------------- * FORMAT STATEMENTS FOR INTERMEDIATE OUTPUT AND ENTRY ERRORS *------------------------------------------------------------------- 1000 FORMAT (//' pass #1 successful:') 1010 FORMAT (' computed pass bounds - hadamard: ',I6) 1020 FORMAT (' modified: ',I6,/) 1025 FORMAT (' hadamard bound could not be satisfied') 1030 FORMAT (//' pass #1 successful, matrix is full rank') 1035 FORMAT (' modified hadamard bound could not be satisfied') 1040 FORMAT (/' pass #1 computed the following',I5, * ' dependent rows ') 1045 FORMAT (/' pass #1 computed the following ',I4, * ' dependent rows in ',/,' the overdetermined system') 1050 FORMAT (' completed pass number ',I5,'. time: ',F12.4) 1060 FORMAT (//' input matrix dimensions are',2I6,/ * ' pass #1:',/) 1080 FORMAT (//' rank has increased with prime =',I12/ * ' discarding all previous primes'//) 1090 FORMAT (//' rank has decreased with prime =',I12/ * ' discarding this prime'//) 1095 FORMAT (//' solution hash table 100% full - cannot', * ' continue passes.',/ * ' continue to output rank and row dependencies') 1100 FORMAT (/5X,'matrix file unit error - unit not opened') 1110 FORMAT (/5X,'output file unit error - unit not opened') 1120 FORMAT (/5X,'intermediate output parameter error - allowed' * ' range (0-4)') 1130 FORMAT (/5X,'intermediate file unit error - unit not opened') 1140 FORMAT (/5X,'output parameter error - allowed range (0-3)') 1150 FORMAT (/5X,'base error - must be less than ',I10) 1160 FORMAT (/5X,'stop criterion error - allowed range > -2') 1170 FORMAT (/5X,'iwrk improperly dimensioned - must be at least ',I6) 1180 FORMAT (/5X, * '********************************************************',/5X, * ' ',/5X, * ' input error(s) encountered in rsdana', * ' ',/5X) 1190 FORMAT (/5X, * '********************************************************',/) 1200 FORMAT (/5X,'hwrk improperly dimensioned - must be at least ',I6) 1205 FORMAT (/5X,'iwrk improperly dimensioned - must be at least ',I6) 1210 FORMAT (/5X,'number of primes too large - max allowed = ',I6) 1215 FORMAT (/5X,'nstor too small - must be at least 2') 1220 FORMAT (/5X,'matrix must have a positive number of rows') 1230 FORMAT (/5X,'matrix must have a positive number of columns') 1240 FORMAT (/5X,'matrix must have a positive number of rhs for', * ' a solution') 1250 FORMAT (/5X,'terminal unit number error ') 1260 FORMAT (/5X,'output style must equal 0 or 1') 1270 FORMAT (/5X,'input file contains multiple ref. to a matrix', * ' element',2I5) 1300 FORMAT (//' processing the particular solutions:',/) 1310 FORMAT (' completed solution: ',I4,'. time: ',F12.4) 1320 FORMAT (//' processing the null space basis, dimension = ',I6) 1330 FORMAT (' completed element: ',I4,'. time: ',F12.4) 1340 FORMAT (/5X, * '********************************************************',/5X, * ' ',/5X, * ' matrix element out of range (i,j) = ',2I6,/5X, * ' ',/5X, * '********************************************************',//) *------------------------------------------------------------------- * FORMAT STATEMENTS FOR HARD COPY RESULTS *------------------------------------------------------------------- 9000 FORMAT (/5X, * '********************************************************',/5X, * ' ',/5X, * ' numerical results ',/5X, * ' ',/5X, * '********************************************************',//) 9010 FORMAT (5X,'input matrix contains: ',I4,' rows',/5X, * ' ',I4,' columns',/5X, * ' ',I4,' rhs',//) 9015 FORMAT (5X,'hadamard bound on passes is: ',I4) 9020 FORMAT (5X,'modified bound on passes is: ',I4) 9025 FORMAT (5X,'passes for termination: ',I4) 9030 FORMAT (5X,'matrix rank is ',I4//,5X, * 'the dependent rows are'/) 9035 FORMAT (5X,'matrix rank is ',I4//,5X, * 'no dependent rows were found'/) 9040 FORMAT (/5X,'the row space is spanned by rows'/) 9045 FORMAT (5X,'unused passes:',I5,//) 9050 FORMAT (/5X,'the column space is spanned by columns'/) 9060 FORMAT (//5X, * '********************************************************',/5X, * ' ',/5X, * ' each of the integers below is in fixed-radix form.',/5X, * ' form left to right, the numbers are:',/5X, * ' i := sign * (coef(1)*base**(l-1) + coef(2)*base**l',/5X, * ' + ... + coef(l-1)*base + coef(l) )',/5X, * ' with base = ',I11,/5X, * ' ',/5X, * '********************************************************',/) 9065 FORMAT (//5X, * '********************************************************',/5X, * ' ',/5X, * ' each of the integers below is in fixed-radix form.',/5X, * ' ',/5X, * '********************************************************',/) 9070 FORMAT (5X,'the minor of the leading non-singular submatrix:',/) 9075 FORMAT (5X,'the determinant of the integer matrix:',/) 9080 FORMAT (45(5X,10I6,/)) 9085 FORMAT (/5X,'the determinant of the rational matrix:',/5X, * 'numerator ',/) 9090 FORMAT (' ') 9100 FORMAT (/5X, * '********************************************************',/5X, * ' ',/5X, * ' non-zero elements in the null space basis ',/5X, * ' ',/5X, * ' read as (i,j) element in pivoted matrix followed ',/5X, * ' the value expressed in fixed radix form ',/5X, * ' ',/5X, * '********************************************************',/) 9110 FORMAT (/5X, * '********************************************************',/5X, * ' ',/5X, * ' null space basis is empty ',/5X, * ' ',/5X, * '********************************************************',/) 9120 FORMAT (/,5X,'null space vector: ',I4,' hash column key: ',I3,/ * 5X,'corresponding matrix column: ',I4/) 9130 FORMAT (/5X,'the null space is ',I6,' dimensional.',/ * 5X,'the corresponding matrix columns are:',/) 9140 FORMAT (3X,'element:',I7) 9200 FORMAT (/5X, * '********************************************************',/5X, * ' ',/5X, * ' solutions to the linear systems ',/5X, * ' ',/5X, * '********************************************************',/) 9210 FORMAT (//5X,'rhs #',I2,' is not in the span of the matrix', * ' column space.',/5X,'thus, no explicit solution exists.') 9220 FORMAT (//5X,'solution to rhs #',I2,'; hash column key ',I3) 9260 FORMAT (/,5X,'numerators ',/) 9270 FORMAT (/,5X,'denominator ',/) 9300 FORMAT (//,5X, * '********************************************************',/5X, * ' ',/5X, * ' routine statistics ',/5X, * ' ',/5X, * ' note: optimal efficiency ratings = 1 ',/5X, * ' ',/5X, * '********************************************************',/) 9310 FORMAT (/,5X, * ' (# tot hashes) + (# double hashes)',/5X, * ' hashing efficiency = ----------------------------------',/5X, * ' (# tot hashes) ',//5X, * ' = ',F12.5) 9320 FORMAT (/,5X, * ' max # of non-zero elements ',/5X, * ' elimination efficiency = ------------------------------',/5X, * ' initial # of non-zero elements',//5X, * ' = ',F12.5) 9330 FORMAT (/5X,' elimination maximum double hashes necessary: ',I6) 9335 FORMAT (/5X,' solution maximum double hashes necessary: ',I6) 9340 FORMAT (/5X,' elimination hash table use: ',I8,' of',I8,' = ' * F6.2,'%') 9345 FORMAT (/5X,' solution hash table use: ',I8,' of',I8,' = ' * F6.2,'%') 9350 FORMAT (/,5X, * ' # initial non-zero + # of fillin elements' * ,/5X, * ' pivot efficiency = -----------------------------------------' * ,/5X, * ' initial # of non-zero elements',//5X, * ' = ',F12.5) 9355 FORMAT (/5X,' output hash table is full - all results are not', * 'stored') 9360 FORMAT (/5X,' output hash table use: ',I8,' of',I8,' = ' * F6.2,'%') 9370 FORMAT (/5X,' suggested minimum for size of work space in', * ' analyz = ',I6) 9380 FORMAT (/5X,' estimate of the required work space in analyz = ', * I7) END