SUBROUTINE MOVE2(NMFILE,NR,JSTART,NJSAVE,JSAVE,T) DIMENSION NMFILE(1),JSAVE(1),T(1) C C THIS SUBROUTINE MOVES THE RECORDS SPECIFIED BY JSAVE C TO THE FRONT OF THE DATA FILE. NJSAVE IS THE C NUMBER OF RECORDS TO BE MOVED AND NMFILE CONTAINS C THE NAME OF THE DATA FILE. C IF(NJSAVE - 1) 40,25,10 C C RANK THE ELEMENTS OF JSAVE IN ASCENDING ORDER. C 10 IJK = 0 DO 20 J=2,NJSAVE J1 = J - 1 IF(JSAVE(J1) .LE. JSAVE(J)) GO TO 20 IJK = JSAVE(J1) JSAVE(J1) = JSAVE(J) JSAVE(J) = IJK 20 CONTINUE IF(IJK .NE. 0) GO TO 10 C C MOVE THE RECORDS SPECIFIED BY JSAVE TO THE BEGINNING C OF THE DATA FILE. C 25 NR1 = NR CALL FDBSET(1,'OLD') CALL ASSIGN(1,NMFILE,10,IER) DEFINE FILE 1(NR1,252,U,IER) DO 30 M=1,NJSAVE MM = JSAVE(M) MMM = M + JSTART - 1 MMMM = MM + JSTART - 1 READ(1'MMMM) (T(I),I=1,126) WRITE(1'MMM) (T(I),I=1,126) 30 CONTINUE CALL CLOSE(1) 40 RETURN END RECORDS SPECIFIED BY JSAVE TO THE BEGINNING C OF THE DATA FILE. C 25 NR1 = NR CALL FDBSET(1,'OLD') CALL ASSIGN(1,NMFILE,10,IER) DEFINE FILE 1(NR1,252,U,IER) DO CV = 100.0 * SD / U1 PRINT 105, SD,CV 105 FORMAT(' STANDARD DEVIATION = ',F14.5,3X,/, 1 ' COEFFICIENT OF VARIANCE = ',F14.5) PRINT 110 110 FORMAT(//,6X,'TAG',12X,'VALUE',7X,'LOG', 1 9X,'RATIO',/) 1 R12 = ALOG(SS1/SS2) C C WRITE THE MEANS AND VARIANCES INTO A SAVE FILE. C CALL FDBSET(1,'APPEND') CALL ASSIGN(1,NMSVFL,10,IER) WRITE(1,111) R12,U1,SS1,U2,SS2 111 FORMAT(1X,5G12.4) CALL CLOSE(1) C C OPEN THE CLASSIFICATION FILE AND READ NJSAVE AND JSAVE C WHICH ARE THE NUMBER AND LIST OF DESCRIPTORS RETAINED C BY THE NONPARAMETRIC SUPERVISED LEARNING ROUTINE. C CALL FDBSET(3,'OLD') CALL ASSIGN(3,NMSVFL,10,IER) READ(3,90) NJSAVE 90 FORMAT(1X,12I5) READ(3,90) (JSAVE(I),I=1,NJSAVE) C C CALL MOVE2 TO MOVE THE DATA CORRESPONDING TO THE LISTED C DESCRIPTORS TO THE BEGINNING OF THE DATA FILE. C CALL MOVE2(NMFILE,NR,JSTART,NJSAVE,JSAVE,T) C C READ FROM THE CLASSIFICATION FILE NJSAVE AND JSAVE WHICH C ARE NOW THE NUMBER AND LIST OF DESCRIPTORS RETAINED C BY THE PARAMETRIC SUPERVISED LEARNING ROUTINE. ALSO READ C THE WEIGHTS W(I) GIVEN TO EACH DESCRIPTOR. C READ(3,90) NJSAVE READ(3,90) (JSAVE(I),I=1,NJSAVE) READ(3,100) (W(I),I=1,NJSAVE) 100 FORMAT(1X,5G12.4) C C CALL MOVE2 TO MOVE THE DATA CORRESPONDING TO THE LISTED C WEIGHTS TO THE BEGINNING OF THE DATA FILE. C CALL MOVE2(NMFILE,NR,JSTART,NJSAVE,JSAVE,T) C C CALL LINCM3 TO CREATE A LINEAR COMBINATION OF DESCRIPTORS Z. C CALL LINCM3(NMFILE,NR,N,JSTART,NJSAVE,W,Z,T) C C READ THE MEANS, U1 AND U2, AND THE VARIANCES, SS1 AND C SS2, FOR THE TWO TRAINING SETS USED IN THE SUPERVISED C LEARNING. ALSO READ R12 WHICH IS DEFINED BY C R12 = ALOG(SS1/SS2) C READ(3,100) R12,U1,SS1,U2,SS2 CALL CLOSE(3) PRINT 110 110 FORMAT(//,4X,'TAG',12X,'VALUE',7X,'LOG', 1 9X,'RATIO',4X,'CLASS') C C FOR EACH CELL EVALUATE THE LOG OF THE LIKELIHOOD RATIO, C TEST; THE LIKELIHOOD RATIO, RVALUE; AND PRINT C THE RESULTS. ALSO CALCULATE THE AVERAGE LOG OF THE C LIKELIHOOD RATIO, AVGL. C AVGL = 0.0 DO 140 I=1,N TEST = RATIO(Z(I)) AVGL = AVGL + TEST IF(TEST .LT. 38.0) GO TO 115 RVALUE = 1.0E38 GO TO 120 115 RVALUE = 10.0 ** TEST 120 II = 1 C C IF TEST IF HEGATIVE CLASSIFY THE CELL IN CLASS 1; IF C IT IS POSITIVE CLASSIFY THE CELL IN CLASS2, C IF(TEST .GT. 0.0) II = 2 PRINT 130,IFNM1,IFNM2,IFNM3,I,Z(I), 1 TEST,RVALUE,II 130 FORMAT(1X,3A2,I5,3X,F10.3,2X,F10.4, 1 2X,E10.2,3X,I1) 140 CONTINUE AVGL = AVGL / N C C DISPLAY THE AVERAGE LOG OF THE LIKELIHOOD RATIO. C PRINT 150, AVGL 150 FORMAT(//,' AVGL = ',F10.3) WRITE(6,160) 160 FORMAT('0END OF RUN FOR CLASFY') CALL SETEF(50,IDS) END 0.0) II = 2 PRINT 130,IFNM1,IFNM2,IFNM3,I,Z(I), 1 TEST,RVALUE,II 130 FORMAT(1X,3A2,I5,3X,F10.3,2X,F10.4, 1 2X,E10.2,3X,I1) 140 CONTINUE AVGL = AVGL / N C C DISPLAY THE AVERAGE LOG OF THE LIKELIHOOD RATIO. C PRINT 150, AVGL 1TIONS HAVE BEEN MADE, RETURN. IF(ITER .GE. ITMAX) RETURN C GET THE GRADIENT OF THE ERROR FUNCTION AT THE FINAL C WEIGHTS. CALL GRAD(N1,N2,JSTART,NCOORD,W, 1 Z1,Z2,T,YERR,DSTEP,WP) C GET THE DIFFERENCE BETWEEN THE INITIAL AND FINAL C GRADIENT. DO 210 I=1,NCOORD 210 YP(I) = WP(I) - WPT(I) C GET A NEW HH MATRIX. CALL MTRX(NCOORD,YP,S) GO TO 20 END EEN MADE, RETURN. IF(ITER .GE. ITMAX) RETURN C GET THE GRADIENT OF THE ERROR FUNCTION AT THE FINAL C WEIGHTS. CALL GRAD(N1,N2,JSTART,NCOORD,W, 1 Z1,Z2,T,YERR,DSTEP,WP) C GET THE DIFFERENCE BETWEEN THE INITIAL AND FINAL C GRADIENT. DO 210 I=1,NCOORD 210 YP(I) = WP(I) - WPT(I) C GET A NEW HH MATRIX. CALL MTRX(NCOORD,YP,S) GO TO 22 (NMSVFL(I),I=1,5) CALL CLOSE(1) JSTART = 1 C READ THE DESCRIPTOR NAMES AND PLACE THEM IN THE ARRAYS C NAME1, NAME2, AND NAME3. CALL FDBSET(1,'OLD') CALL ASSIGN(1,NMFIL3,10,IER) DEFINE FILE 1(NR1,3,U,IER) DO 60 I=1,NCORD1 60 READ(1'I) NAME1(I),NAME2(I),NAME3(I) CALL CLOSE(1) C C COPY THE CELL DATA INTO THE FILES VECTR1.DAT AND C VECTR2.DAT. C CALL FDBSET(1,'OLD') CALL ASSIGN(1,NMFIL1,10,IER) DEFINE FILE 1(NR1,252,U,IER) CALL FDBSET(2,'NEW') CALL ASSIGN(2,'VECTR1.DAT',10,IER) DEFINE FILE 2(24,252,U,IER) DO 65 I=1,NCORD1 READ(1'I) (T(J),J=1,126) I1 = JSTART + I - 1 WRITE(2'I1) (T(J),J=1,126) 65 CONTINUE CALL CLOSE(1) CALL CLOSE(2) C CALL FDBSET(1,'OLD') CALL ASSIGN(1,NMFIL2,10,IER) DEFINE FILE 1(NR1,252,U,IER) CALL FDBSET(2,'NEW') CALL ASSIGN(2,'VECTR2.DAT',10,IER) DEFINE FILE 2(24,252,U,IER) DO 70 I=1,NCORD1 READ(1'I) (T(J),J=1,126) I1 = JSTART + I - 1 WRITE(2'I1) (T(J),J=1,126) 70 CONTINUE CALL CLOSE(1) CALL CLOSE(2) C C CALL EXTRT1 TO GET THE WEIGHTS W(I) FOR EACH DESCRIPTOR C AND THE LINEAR COMBINATIONS Z1 AND Z2. C CALL EXTRT1(JSTART,NCORD1,N1,N2,W,PW,WT,JSAVE, 2 Z1,Z2,T) C PRINT THE CHOSEN DESCRIPTORS AND THEIR CORRESPONDING C WEIGHTS. PRINT 80,NCORD1 80 FORMAT(////,' THE BEST ',I2,' DESCRIPTORS CHOSEN BY ', 1 'EXTRACT ARE:',//,' DESCRIPTOR',4X,'NAME', 2 7X,'WEIGHT',5X,'RELATIVE WGT',//) DO 85 I=1,NCORD1 J = JSAVE(I) 85 PRINT 90,J,NAME1(J),NAME2(J),NAME3(J),W(I),PW(I) 90 FORMAT(4X,I3,7X,3A2,3X,G12.4,4X,F7.4,//) C C WRITE INTO THE SAVE FILE THE DESCRIPTORS SAVED C AND THEIR CORRESPONDING WEIGHTS. CALL FDBSET(1,'APPEND') CALL ASSIGN(1,NMSVFL,10,IER) WRITE(1,92) NCORD1 92 FORMAT(1X,12I5) WRITE(1,92) (JSAVE(I),I=1,NCORD1) WRITE(1,95) (W(I),I=1,NCORD1) 95 FORMAT(1X,5G12.4) CALL CLOSE(1) C C WRITE THE NEW DESCRIPTOR DATA, Z1(I) AND Z2(I), INTO C THE CELL DATA FILES, VECTR1.DAT AND VECTR2.DAT. CALL FDBSET(1,'OLD') CALL ASSIGN(1,'VECTR1.DAT',10,IER) DEFINE FILE 1(24,252,U,IER) WRITE(1'JSTART) (Z1(I),I=1,126) CALL CLOSE(1) CALL FDBSET(1,'OLD') CALL ASSIGN(1,'VECTR2.DAT',10,IER) DEFINE FILE 1(24,252,U,IER) WRITE(1'JSTART) (Z2(I),I=1,126) CALL CLOSE(1) WRITE(6,100) 100 FORMAT('0END OF FUN FOR STEP 2(EXTRACT)') CALL SETEF(50,I) END LE 1(24,252,U,IER) WRITE(1'JSTART) (Z1(I),I=1,126) DSUM=0.D0 IDO=M DO 41 I=1,M HELP=PIV(IDO) IF(K) 39,39,40 39 HELP=-HELP 40 DSUM=DSUM+DBLE(HELP) PIV(IDO+1)=HELP 41 IDO=IDO-1 PIV(L)=-DSUM PIV(1)=1. C C TRANSFORM VECTOR PIV WITH ROWS OF MATRIX T IDO=IND IF(ILAB) 44,44,42 42 K=1 43 IDO=K 44 DSUM=0.D0 HELP=0. C C START MULTIPLICATION-LOOP DO 46 I=1,L DSUM=DSUM+DBLE(PIV(I)*T(IDO)) TOL=ABS(SNGL(DSUM)) IF(TOL-HELP) 46,46,45 45 HELP=TOL 46 IDO=IDO+L C END OF MULTIPLICATION-LOOP C TOL=1.E-5*HELP IF(ABS(SNGL(DSUM))-TOL) 47,47,48 47 DSUM=0.D0 48 IF(ILAB) 51,51,49 49 I=K+L PIV(I)=DSUM C C TEST FOR LAST COLUMN-TERM K=K+1 IF(K-L) 43,43,18 50 I=(K-1)*L+IND DSUM=T(I) C C COMPUTE NEW TOP-ELEMENT 51 DSUM=DSUM*DBLE(SAVE) TOL=1.E-5*ABS(SNGL(DSUM)) TOP(J)=TOP(J)+SNGL(DSUM) IF(ABS(TOP(J))-TOL) 52,52,53 52 TOP(J)=0. C C TEST FOR LAST TOP-TERM 53 J=J-1 IF(J) 54,54,33 C END OF TOP-LOOP C C TRANSFORM PIVOT-COLUMN 54 I=IND+L PIV(I)=-1. DO 55 I=1,L J=I+L 55 PIV(I)=-PIV(J)*REPI C C UPDATE TRANSFORMATION-MATRIX T J=0 DO 57 I=1,L IDO=J+IND SAVE=T(IDO) T(IDO)=0. DO 56 K=1,L ISE=K+J 56 T(ISE)=T(ISE)+SAVE*PIV(K) 57 J=J+L C C UPDATE INDEX-VECTOR IHE C INITIALIZE CHARACTERISTICS J=0 K=0 ISE=0 IDO=0 C C START QUESTION-LOOP DO 61 I=1,L LL=I+L ILAB=IHE(LL) IF(IHE(I)-IPIV) 59,58,59 58 ISE=I J=ILAB 59 IF(ILAB-IND) 61,60,61 60 IDO=I K=IHE(I) 61 CONTINUE C END OF QUESTION-LOOP C C START MODIFICATION IF(K) 62,62,63 62 IHE(IDO)=IPIV IF(ISE) 67,67,65 63 IF(IND-J) 64,66,64 64 LL=L+L+L+NAN K=K+LL I=IPIV+LL ILAB=IHE(K) IHE(K)=IHE(I) IHE(I)=ILAB IF(ISE) 67,67,65 65 IDO=IDO+L I=ISE+L IHE(IDO)=J IHE(I)=IND 66 IHE(ISE)=0 67 LL=L+L J=LL+IND I=LL+L+IPIV ILAB=IHE(I) IHE(I)=IHE(J) IHE(J)=ILAB C END OF MODIFICATION C GO TO 8 C C SET ERROR PARAMETER IER=-1 SINCE NO SUITABLE PIVOT IS FOUND 68 IER=-1 C C EVALUATE FINAL TABLEAU C COMPUTE SAVE AS MAXIMUM ERROR OF APPROXIMATION AND C HELP AS ADDITIVE CONSTANCE FOR RESULTING COEFFICIENTS 69 SAVE=0. HELP=0. K=L+L+L DO 73 I=1,NAN IDO=K+I J=IHE(IDO) IF(J) 71,70,73 70 SAVE=-TOP(I) 71 IF(M+J+1) 73,72,73 72 HELP=TOP(I) 73 CONTINUE C C PREPARE T,TOP,PIV T(1)=SAVE IDO=NAN+1 J=NAN+N DO 74 I=IDO,J 74 TOP(I)=SAVE DO 75 I=1,M 75 PIV(I)=HELP C C COMPUTE COEFFICIENTS OF RESULTING POLYNOMIAL IN PIV(1) UP TO PI C AND CALCULATE ERRORS AT GIVEN NODES IN TOP(1) UP TO TOP(N) DO 79 I=1,NAN IDO=K+I J=IHE(IDO) IF(J) 76,79,77 76 J=-J PIV(J)=HELP-TOP(I) GO TO 79 77 IF(J-N) 78,78,79 78 J=J+NAN TOP(J)=SAVE+TOP(I) 79 CONTINUE DO 80 I=1,N IDO=NAN+I 80 TOP(I)=TOP(IDO) 81 RETURN END C C .................................................................. C C SUBROUTINE ARAT C C PURPOSE C CALCULATE BEST RATIONAL APPROXIMATION OF A DISCRETE C FUNCTION IN THE LEAST SQUARES SENSE C C USAGE C CALL ARAT(DATI,N,WORK,P,IP,IQ,IER) C C DESCRIPTION OF PARAMETERS C DATI - TWODIMENSIONAL ARRAY WITH 3 COLUMNS AND N ROWS C THE FIRST COLUMN MUST CONTAIN THE GIVEN ARGUMENTS, C THE SECOND COLUMN THE GIVEN FUNCTION VALUES AND C THE THIRD COLUMN THE GIVEN WEIGHTS IF ANY. C IF NO WEIGHTS ARE TO BE USED THEN THE THIRD C COLUMN MAY BE DROPPED , EXCEPT THE FIRST ELEMENT C WHICH MUST CONTAIN A NONPOSITIVE VALUE C N - NUMBER OF NODES OF THE GIVEN DISCRETE FUNCTION C WORK - WORKING STORAGE WHICH IS OF DIMENSION C (IP+IQ)*(IP+IQ+1)+4*N+1 AT LEAST. C ON RETURN THE VALUES OF THE NUMERATOR ARE CONTAINED C IN WORK(N+1) UP TO WORK(2*N), WHILE THE VALUES OF C THE DENOMINATOR ARE STORED IN WORK(2*N+1) UP TO C WORK(3*N) C P - RESULTANT COEFFICIENT VECTOR OF DENOMINATOR AND C NUMERATOR. THE DENOMINATOR IS STORED IN FIRST IQ C LOCATIONS, THE NUMERATOR IN THE FOLLOWING IP C LOCATIONS. C COEFFICIENTS ARE ORDERED FROM LOW TO HIGH. C IP - DIMENSION OF THE NUMERATOR (INPUT VALUE) C IQ - DIMENSION OF THE DENOMINATOR (INPUT VALUE) C IER - RESULTANT ERROR PARAMETER C IER =-1 MEANS FORMAL ERRORS C IER = 0 MEANS NO ERRORS C IER = 1,2 MEANS POOR CONVERGENCE OF ITERATION C IER IS ALSO USED AS INPUT VALUE C A NONZERO INPUT VALUE INDICATES AVAILABILITY OF AN C INITIAL APPROXIMATION STORED IN P C C REMARKS C THE COEFFICIENT VECTORS OF THE DENOMINATOR AND NUMERATOR C OF THE RATIONAL APPROXIMATION ARE BOTH STORED IN P C STARTING WITH LOW POWERS (DENOMINATOR FIRST). C IP+IQ MUST NOT EXCEED N, ALL THREE VALUES MUST BE POSITIVE. C SINCE CHEBYSHEV POLYNOMIALS ARE USED AS FUNDAMENTAL C FUNCTIONS, THE ARGUMENTS SHOULD BE REDUCED TO THE INTERVAL C (-1,1). THIS CAN ALWAYS BE ACCOMPLISHED BY MEANS OF A LINEAR C TRANSFORMATION OF THE ORIGINALLY GIVEN ARGUMENTS. C IF A FIT IN OTHER FUNCTIONS IS REQUIRED, CNP AND CNPS MUST C BE REPLACED BY SUBROUTINES WHICH ARE OF ANALOGOUS DESIGN. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C APLL, APFS, FRAT, CNPS, CNP C CNP IS REQUIRED WITHIN FRAT C C METHOD C THE ITERATIVE SCHEME USED FOR CALCULATION OF THE C APPROXIMATION IS REPEATED SOLUTION OF THE NORMAL EQUATIONS C WHICH ARE OBTAINED BY LINEARIZATION. C A REFINED TECHNIQUE OF THIS LINEAR LEAST SQUARES APPROACH C IS USED WHICH GUARANTEES THAT THE DENOMINATOR IS FREE OF C ZEROES WITHIN THE APPROXIMATION INTERVAL. C FOR REFERENCE SEE C D.BRAESS, UEBER DAEMPFUNG BEI MINIMALISIERUNGSVERFAHREN, C COMPUTING(1966), VOL.1, ED.3, PP.264-272. C D.W.MARQUARDT, AN ALGORITHM FOR LEAST-SQUARES ESTIMATION C OF NONLINEAR PARAMETERS, C JSIAM(1963), VOL.11, ED.2, PP.431-441. C C .................................................................. C SUBROUTINE ARAT(DATI,N,WORK,P,IP,IQ,IER) C C EXTERNAL FRAT C C DIMENSIONED LOCAL VARIABLE DIMENSION IERV(3) C C DIMENSIONED DUMMY VARIABLES DIMENSION DATI(1),WORK(1),P(1) C C INITIALIZE TESTVALUES LIMIT=20 ETA =1.E-11 EPS=1.E-5 C C CHECK FOR FORMAL ERRORS IF(N)4,4,1 1 IF(IP)4,4,2 2 IF(IQ)4,4,3 3 IPQ=IP+IQ IF(N-IPQ)4,5,5 C C ERROR RETURN IN CASE OF FORMAL ERRORS 4 IER=-1 RETURN C C INITIALIZE ITERATION PROCESS 5 KOUNT=0 IERV(2)=IP IERV(3)=IQ NDP=N+N+1 NNE=NDP+NDP IX=IPQ-1 IQP1=IQ+1 IRHS=NNE+IPQ*IX/2 IEND=IRHS+IX C C TEST FOR AVAILABILITY OF AN INITIAL APPROXIMATION IF(IER)8,6,8 C C INITIALIZE NUMERATOR AND DENOMINATOR 6 DO 7 I=2,IPQ 7 P(I)=0. P(1)=1. C C CALCULATE VALUES OF NUMERATOR AND DENOMINATOR FOR INITIAL C APPROXIMATION 8 DO 9 J=1,N T=DATI(J) I=J+N CALL CNPS(WORK(I),T,P(IQP1),IP) K=I+N 9 CALL CNPS(WORK(K),T,P,IQ) C C SET UP NORMAL EQUATIONS (MAIN LOOP OF ITERATION) 10 CALL APLL(FRAT,N,IX,WORK,WORK(IEND+1),DATI,IERV) C C CHECK FOR ZERO DENOMINATOR IF(IERV(1))4,11,4 11 INCR=0 RELAX=2. C C RESTORE MATRIX IN WORKING STORAGE 12 J=IEND DO 13 I=NNE,IEND J=J+1 13 WORK(I)=WORK(J) IF(KOUNT)14,14,15 C C SAVE SQUARE SUM OF ERRORS 14 OSUM=WORK(IEND) DIAG=OSUM*EPS K=IQ C C ADD CONSTANT TO DIAGONAL IF(WORK(NNE))17,17,19 15 IF(INCR)19,19,16 16 K=IPQ 17 J=NNE-1 DO 18 I=1,K WORK(J)=WORK(J)+DIAG 18 J=J+I C C SOLVE NORMAL EQUATIONS 19 CALL APFS(WORK(NNE),IX,IRES,1,EPS,ETA,IER) C C CHECK FOR FAILURE OF EQUATION SOLVER IF(IRES)4,4,20 C C TEST FOR DEFECTIVE NORMALEQUATIONS 20 IF(IRES-IX)21,24,24 21 IF(INCR)22,22,23 22 DIAG=DIAG*0.125 23 DIAG=DIAG+DIAG INCR=INCR+1 C C START WITH OVER RELAXATION RELAX=8. IF(INCR-LIMIT)12,45,45 C C CALCULATE VALUES OF CHANGE OF NUMERATOR AND DENOMINATOR 24 L=NDP J=NNE+IRES*(IRES-1)/2-1 K=J+IQ WORK(J)=0. IRQ=IQ IRP=IRES-IQ+1 IF(IRP)25,26,26 25 IRQ=IRES+1 26 DO 29 I=1,N T=DATI(I) WORK(I)=0. CALL CNPS(WORK(I),T,WORK(K),IRP) M=L+N CALL CNPS(WORK(M),T,WORK(J),IRQ) IF(WORK(M)*WORK(L))27,29,29 27 SUM=WORK(L)/WORK(M) IF(RELAX+SUM)29,29,28 28 RELAX=-SUM 29 L=L+1 C C MODIFY RELAXATION FACTOR IF NECESSARY SSOE=OSUM ITER=LIMIT 30 SUM=0. RELAX=RELAX*0.5 DO 32 I=1,N M=I+N K=M+N L=K+N SAVE=DATI(M)-(WORK(M)+RELAX*WORK(I))/(WORK(K)+RELAX*WORK(L)) SAVE=SAVE*SAVE IF(DATI(NDP))32,32,31 31 SAVE=SAVE*DATI(K) 32 SUM=SUM+SAVE IF(ITER)45,33,33 33 ITER=ITER-1 IF(SUM-OSUM)34,37,35 34 OSUM=SUM GOTO 30 C C TEST FOR IMPROVEMENT 35 IF(OSUM-SSOE)36,30,30 36 RELAX=RELAX+RELAX 37 T=0. SAVE=0. K=IRES+1 DO 38 I=2,K J=J+1 T=T+ABS(P(I)) P(I)=P(I)+RELAX*WORK(J) 38 SAVE=SAVE+ABS(P(I)) C C UPDATE CURRENT VALUES OF NUMERATOR AND DENOMINATOR DO 39 I=1,N J=I+N K=J+N L=K+N WORK(J)=WORK(J)+RELAX*WORK(I) 39 WORK(K)=WORK(K)+RELAX*WORK(L) C C TEST FOR CONVERGENCE IF(INCR)40,40,42 40 IF(SSOE-OSUM-RELAX*EPS*OSUM)46,46,41 41 IF(ABS(T-SAVE)-RELAX*EPS*SAVE)46,46,42 42 IF(OSUM-ETA*SAVE)46,46,43 43 KOUNT=KOUNT+1 IF(KOUNT-LIMIT)10,44,44 C C ERROR RETURN IN CASE OF POOR CONVERGENCE 44 IER=2 RETURN 45 IER=1 RETURN C C NORMAL RETURN 46 IER=0 RETURN END C C .................................................................. C C SUBROUTINE ARRAY C C PURPOSE C CONVERT DATA ARRAY FROM SINGLE TO DOUBLE DIMENSION OR VICE C VERSA. THIS SUBROUTINE IS USED TO LINK THE USER PROGRAM C WHICH HAS DOUBLE DIMENSION ARRAYS AND THE SSP SUBROUTINES C WHICH OPERATE ON ARRAYS OF DATA IN A VECTOR FASHION. C C USAGE C CALL ARRAY (MODE,I,J,N,M,S,D) C C DESCRIPTION OF PARAMETERS C MODE - CODE INDICATING TYPE OF CONVERSION C 1 - FROM SINGLE TO DOUBLE DIMENSION C 2 - FROM DOUBLE TO SINGLE DIMENSION C I - NUMBER OF ROWS IN ACTUAL DATA MATRIX C J - NUMBER OF COLUMNS IN ACTUAL DATA MATRIX C N - NUMBER OF ROWS SPECIFIED FOR THE MATRIX D IN C DIMENSION STATEMENT C M - NUMBER OF COLUMNS SPECIFIED FOR THE MATRIX D IN C DIMENSION STATEMENT C S - IF MODE=1, THIS VECTOR IS INPUT WHICH CONTAINS THE C ELEMENTS OF A DATA MATRIX OF SIZE I BY J. COLUMN I+1 C OF DATA MATRIX FOLLOWS COLUMN I, ETC. IF MODE=2, C THIS VECTOR IS OUTPUT REPRESENTING A DATA MATRIX OF C SIZE I BY J CONTAINING ITS COLUMNS CONSECUTIVELY. C THE LENGTH OF S IS IJ, WHERE IJ=I*J. C D - IF MODE=1, THIS MATRIX OF SIZE N BY M IS OUTPUT, C CONTAINING A DATA MATRIX OF SIZE I BY J IN THE FIRST C I ROWS AND J COLUMNS. IF MODE=2, THIS N BY M MATRIX C IS INPUT CONTAINING A DATA MATRIX OF SIZE I BY J IN C THE FIRST I ROWS AND J COLUMNS. C C REMARKS C VECTOR S CAN BE IN THE SAME LOCATION AS MATRIX D. VECTOR S C IS REFERRED AS A MATRIX IN OTHER SSP ROUTINES, SINCE IT C CONTAINS A DATA MATRIX. C THIS SUBROUTINE CONVERTS ONLY GENERAL DATA MATRICES (STORAGE C MODE OF 0). C C SUBROUTINES AND FUNCTION SUBROUTINES REQUIRED C NONE C C METHOD C REFER TO THE DISCUSSION ON VARIABLE DATA SIZE IN THE SECTION C DESCRIBING OVERALL RULES FOR USAGE IN THIS MANUAL. C C .................................................................. C SUBROUTINE ARRAY (MODE,I,J,N,M,S,D) DIMENSION S(1),D(1) C NI=N-I C C TEST TYPE OF CONVERSION C IF(MODE-1) 100, 100, 120 C C CONVERT FROM SINGLE TO DOUBLE DIMENSION C 100 IJ=I*J+1 NM=N*J+1 DO 110 K=1,J NM=NM-NI DO 110 L=1,I IJ=IJ-1 NM=NM-1 110 D(NM)=S(IJ) GO TO 140 C C CONVERT FROM DOUBLE TO SINGLE DIMENSION C 120 IJ=0 NM=0 DO 130 K=1,J DO 125 L=1,I IJ=IJ+1 NM=NM+1 125 S(IJ)=D(NM) 130 NM=NM+NI C 140 RETURN END C C .................................................................. C C SUBROUTINE ATEIG C C PURPOSE C COMPUTE THE EIGENVALUES OF A REAL ALMOST TRIANGULAR MATRIX C C USAGE C CALL ATEIG(M,A,RR,RI,IANA,IA) C C DESCRIPTION OF THE PARAMETERS C M ORDER OF THE MATRIX C A THE INPUT MATRIX, M BY M C RR VECTOR CONTAINING THE REAL PARTS OF THE EIGENVALUES C ON RETURN C RI VECTOR CONTAINING THE IMAGINARY PARTS OF THE EIGEN- C VALUES ON RETURN C IANA VECTOR WHOSE DIMENSION MUST BE GREATER THAN OR EQUAL C TO M, CONTAINING ON RETURN INDICATIONS ABOUT THE WAY C THE EIGENVALUES APPEARED (SEE MATH. DESCRIPTION) C IA SIZE OF THE FIRST DIMENSION ASSIGNED TO THE ARRAY A C IN THE CALLING PROGRAM WHEN THE MATRIX IS IN DOUBLE C SUBSCRIPTED DATA STORAGE MODE. C IA=M WHEN THE MATRIX IS IN SSP VECTOR STORAGE MODE. C C REMARKS C THE ORIGINAL MATRIX IS DESTROYED C THE DIMENSION OF RR AND RI MUST BE GREATER OR EQUAL TO M C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C QR DOUBLE ITERATION C C REFERENCES C J.G.F. FRANCIS - THE QR TRANSFORMATION---THE COMPUTER C JOURNAL, VOL. 4, NO. 3, OCTOBER 1961, VOL. 4, NO. 4, JANUARY C 1962. J. H. WILKINSON - THE ALGEBRAIC EIGENVALUE PROBLEM - C CLARENDON PRESS, OXFORD, 1965. C C .................................................................. C SUBROUTINE ATEIG(M,A,RR,RI,IANA,IA) DIMENSION A(1),RR(1),RI(1),PRR(2),PRI(2),IANA(1) INTEGER P,P1,Q C E7=1.0E-8 E6=1.0E-6 E10=1.0E-10 DELTA=0.5 MAXIT=30 C C INITIALIZATION C N=M 20 N1=N-1 IN=N1*IA NN=IN+N IF(N1) 30,1300,30 30 NP=N+1 C C ITERATION COUNTER C IT=0 C C ROOTS OF THE 2ND ORDER MAIN SUBMATRIX AT THE PREVIOUS C ITERATION C DO 40 I=1,2 PRR(I)=0.0 40 PRI(I)=0.0 C C LAST TWO SUBDIAGONAL ELEMENTS AT THE PREVIOUS ITERATION C PAN=0.0 PAN1=0.0 C C ORIGIN SHIFT C R=0.0 S=0.0 C C ROOTS OF THE LOWER MAIN 2 BY 2 SUBMATRIX C N2=N1-1 IN1=IN-IA NN1=IN1+N N1N=IN+N1 N1N1=IN1+N1 60 T=A(N1N1)-A(NN) U=T*T V=4.0*A(N1N)*A(NN1) IF(ABS(V)-U*E7) 100,100,65 65 T=U+V IF(ABS(T)-AMAX1(U,ABS(V))*E6) 67,67,68 67 T=0.0 68 U=(A(N1N1)+A(NN))/2.0 V=SQRT(ABS(T))/2.0 IF(T)140,70,70 70 IF(U) 80,75,75 75 RR(N1)=U+V RR(N)=U-V GO TO 130 80 RR(N1)=U-V RR(N)=U+V GO TO 130 100 IF(T)120,110,110 110 RR(N1)=A(N1N1) RR(N)=A(NN) GO TO 130 120 RR(N1)=A(NN) RR(N)=A(N1N1) 130 RI(N)=0.0 RI(N1)=0.0 GO TO 160 140 RR(N1)=U RR(N)=U RI(N1)=V RI(N)=-V 160 IF(N2)1280,1280,180 C C TESTS OF CONVERGENCE C 180 N1N2=N1N1-IA RMOD=RR(N1)*RR(N1)+RI(N1)*RI(N1) EPS=E10*SQRT(RMOD) IF(ABS(A(N1N2))-EPS)1280,1280,240 240 IF(ABS(A(NN1))-E10*ABS(A(NN))) 1300,1300,250 250 IF(ABS(PAN1-A(N1N2))-ABS(A(N1N2))*E6) 1240,1240,260 260 IF(ABS(PAN-A(NN1))-ABS(A(NN1))*E6)1240,1240,300 300 IF(IT-MAXIT) 320,1240,1240 C C COMPUTE THE SHIFT C 320 J=1 DO 360 I=1,2 K=NP-I IF(ABS(RR(K)-PRR(I))+ABS(RI(K)-PRI(I))-DELTA*(ABS(RR(K)) 1 +ABS(RI(K)))) 340,360,360 340 J=J+I 360 CONTINUE GO TO (440,460,460,480),J 440 R=0.0 S=0.0 GO TO 500 460 J=N+2-J R=RR(J)*RR(J) S=RR(J)+RR(J) GO TO 500 480 R=RR(N)*RR(N1)-RI(N)*RI(N1) S=RR(N)+RR(N1) C C SAVE THE LAST TWO SUBDIAGONAL TERMS AND THE ROOTS OF THE C SUBMATRIX BEFORE ITERATION C 500 PAN=A(NN1) PAN1=A(N1N2) DO 520 I=1,2 K=NP-I PRR(I)=RR(K) 520 PRI(I)=RI(K) C C SEARCH FOR A PARTITION OF THE MATRIX, DEFINED BY P AND Q C P=N2 IF (N-3)600,600,525 525 IPI=N1N2 DO 580 J=2,N2 IPI=IPI-IA-1 IF(ABS(A(IPI))-EPS) 600,600,530 530 IPIP=IPI+IA IPIP2=IPIP+IA D=A(IPIP)*(A(IPIP)-S)+A(IPIP2)*A(IPIP+1)+R IF(D)540,560,540 540 IF(ABS(A(IPI)*A(IPIP+1))*(ABS(A(IPIP)+A(IPIP2+1)-S)+ABS(A(IPIP2+2) 1 )) -ABS(D)*EPS) 620,620,560 560 P=N1-J 580 CONTINUE 600 Q=P GO TO 680 620 P1=P-1 Q=P1 IF (P1-1) 680,680,650 650 DO 660 I=2, P1 IPI=IPI-IA-1 IF(ABS(A(IPI))-EPS)680,680,660 660 Q=Q-1 C C QR DOUBLE ITERATION C 680 II=(P-1)*IA+P DO 1220 I=P,N1 II1=II-IA IIP=II+IA IF(I-P)720,700,720 700 IPI=II+1 IPIP=IIP+1 C C INITIALIZATION OF THE TRANSFORMATION C G1=A(II)*(A(II)-S)+A(IIP)*A(IPI)+R G2=A(IPI)*(A(IPIP)+A(II)-S) G3=A(IPI)*A(IPIP+1) A(IPI+1)=0.0 GO TO 780 720 G1=A(II1) G2=A(II1+1) IF(I-N2)740,740,760 740 G3=A(II1+2) GO TO 780 760 G3=0.0 780 CAP=SQRT(G1*G1+G2*G2+G3*G3) IF(CAP)800,860,800 800 IF(G1)820,840,840 820 CAP=-CAP 840 T=G1+CAP PSI1=G2/T PSI2=G3/T ALPHA=2.0/(1.0+PSI1*PSI1+PSI2*PSI2) GO TO 880 860 ALPHA=2.0 PSI1=0.0 PSI2=0.0 880 IF(I-Q)900,960,900 900 IF(I-P)920,940,920 920 A(II1)=-CAP GO TO 960 940 A(II1)=-A(II1) C C ROW OPERATION C 960 IJ=II DO 1040 J=I,N T=PSI1*A(IJ+1) IF(I-N1)980,1000,1000 980 IP2J=IJ+2 T=T+PSI2*A(IP2J) 1000 ETA=ALPHA*(T+A(IJ)) A(IJ)=A(IJ)-ETA A(IJ+1)=A(IJ+1)-PSI1*ETA IF(I-N1)1020,1040,1040 1020 A(IP2J)=A(IP2J)-PSI2*ETA 1040 IJ=IJ+IA C C COLUMN OPERATION C IF(I-N1)1080,1060,1060 1060 K=N GO TO 1100 1080 K=I+2 1100 IP=IIP-I DO 1180 J=Q,K JIP=IP+J JI=JIP-IA T=PSI1*A(JIP) IF(I-N1)1120,1140,1140 1120 JIP2=JIP+IA T=T+PSI2*A(JIP2) 1140 ETA=ALPHA*(T+A(JI)) A(JI)=A(JI)-ETA A(JIP)=A(JIP)-ETA*PSI1 IF(I-N1)1160,1180,1180 1160 A(JIP2)=A(JIP2)-ETA*PSI2 1180 CONTINUE IF(I-N2)1200,1220,1220 1200 JI=II+3 JIP=JI+IA JIP2=JIP+IA ETA=ALPHA*PSI2*A(JIP2) A(JI)=-ETA A(JIP)=-ETA*PSI1 A(JIP2)=A(JIP2)-ETA*PSI2 1220 II=IIP+1 IT=IT+1 GO TO 60 C C END OF ITERATION C 1240 IF(ABS(A(NN1))-ABS(A(N1N2))) 1300,1280,1280 C C TWO EIGENVALUES HAVE BEEN FOUND C 1280 IANA(N)=0 IANA(N1)=2 N=N2 IF(N2)1400,1400,20 C C ONE EIGENVALUE HAS BEEN FOUND C 1300 RR(N)=A(NN) RI(N)=0.0 IANA(N)=1 IF(N1)1400,1400,1320 1320 N=N1 GO TO 20 1400 RETURN END C C .................................................................. C C SUBROUTINE ATSE C C PURPOSE C NDIM POINTS OF A GIVEN TABLE WITH EQUIDISTANT ARGUMENTS ARE C SELECTED AND ORDERED SUCH THAT C ABS(ARG(I)-X).GE.ABS(ARG(J)-X) IF I.GT.J. C C USAGE C CALL ATSE (X,ZS,DZ,F,IROW,ICOL,ARG,VAL,NDIM) C C DESCRIPTION OF PARAMETERS C X - THE SEARCH ARGUMENT. C ZS - THE STARTING VALUE OF ARGUMENTS. C DZ - THE INCREMENT OF ARGUMENT VALUES. C F - IN CASE ICOL=1, F IS THE VECTOR OF FUNCTION VALUES C (DIMENSION IROW). C IN CASE ICOL=2, F IS AN IROW BY 2 MATRIX. THE FIRST C COLUMN SPECIFIES THE VECTOR OF FUNCTION VALUES AND C THE SECOND THE VECTOR OF DERIVATIVES. C IROW - THE DIMENSION OF EACH COLUMN IN MATRIX F. C ICOL - THE NUMBER OF COLUMNS IN F (I.E. 1 OR 2). C ARG - THE RESULTING VECTOR OF SELECTED AND ORDERED C ARGUMENT VALUES (DIMENSION NDIM). C VAL - THE RESULTING VECTOR OF SELECTED FUNCTION VALUES C (DIMENSION NDIM) IN CASE ICOL=1. IN CASE ICOL=2, C VAL IS THE VECTOR OF FUNCTION AND DERIVATIVE VALUES C (DIMENSION 2*NDIM) WHICH ARE STORED IN PAIRS (I.E. C EACH FUNCTION VALUE IS FOLLOWED BY ITS DERIVATIVE C VALUE). C NDIM - THE NUMBER OF POINTS WHICH MUST BE SELECTED OUT OF C THE GIVEN TABLE. C C REMARKS C NO ACTION IN CASE IROW LESS THAN 1. C IF INPUT VALUE NDIM IS GREATER THAN IROW, THE PROGRAM C SELECTS ONLY A MAXIMUM TABLE OF IROW POINTS. THEREFORE THE C USER OUGHT TO CHECK CORRESPONDENCE BETWEEN TABLE (ARG,VAL) C AND ITS DIMENSION BY COMPARISON OF NDIM AND IROW, IN ORDER C TO GET CORRECT RESULTS IN FURTHER WORK WITH TABLE (ARG,VAL). C THIS TEST MAY BE DONE BEFORE OR AFTER CALLING C SUBROUTINE ATSE. C SUBROUTINE ATSE ESPECIALLY CAN BE USED FOR GENERATING THE C TABLE (ARG,VAL) NEEDED IN SUBROUTINES ALI, AHI, AND ACFI. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C SELECTION IS DONE BY COMPUTING THE SUBSCRIPT J OF THAT C ARGUMENT, WHICH IS NEXT TO X. C AFTERWARDS NEIGHBOURING ARGUMENT VALUES ARE TESTED AND C SELECTED IN THE ABOVE SENSE. C C .................................................................. C SUBROUTINE ATSE(X,ZS,DZ,F,IROW,ICOL,ARG,VAL,NDIM) C C DIMENSION F(1),ARG(1),VAL(1) IF(IROW-1)19,17,1 C C CASE DZ=0 IS CHECKED OUT 1 IF(DZ)2,17,2 2 N=NDIM C C IF N IS GREATER THAN IROW, N IS SET EQUAL TO IROW. IF(N-IROW)4,4,3 3 N=IROW C C COMPUTATION OF STARTING SUBSCRIPT J. 4 J=(X-ZS)/DZ+1.5 IF(J)5,5,6 5 J=1 6 IF(J-IROW)8,8,7 7 J=IROW C C GENERATION OF TABLE ARG,VAL IN CASE DZ.NE.0. 8 II=J JL=0 JR=0 DO 16 I=1,N ARG(I)=ZS+FLOAT(II-1)*DZ IF(ICOL-2)9,10,10 9 VAL(I)=F(II) GOTO 11 10 VAL(2*I-1)=F(II) III=II+IROW VAL(2*I)=F(III) 11 IF(J+JR-IROW)12,15,12 12 IF(J-JL-1)13,14,13 13 IF((ARG(I)-X)*DZ)14,15,15 14 JR=JR+1 II=J+JR GOTO 16 15 JL=JL+1 II=J-JL 16 CONTINUE RETURN C C CASE DZ=0 17 ARG(1)=ZS VAL(1)=F(1) IF(ICOL-2)19,19,18 18 VAL(2)=F(2) 19 RETURN END C C .................................................................. C C SUBROUTINE ATSG C C PURPOSE C NDIM POINTS OF A GIVEN GENERAL TABLE ARE SELECTED AND C ORDERED SUCH THAT ABS(ARG(I)-X).GE.ABS(ARG(J)-X) IF I.GT.J. C C USAGE C CALL ATSG (X,Z,F,WORK,IROW,ICOL,ARG,VAL,NDIM) C C DESCRIPTION OF PARAMETERS C X - THE SEARCH ARGUMENT. C Z - THE VECTOR OF ARGUMENT VALUES (DIMENSION IROW). C F - IN CASE ICOL=1, F IS THE VECTOR OF FUNCTION VALUES C (DIMENSION IROW). C IN CASE ICOL=2, F IS AN IROW BY 2 MATRIX. THE FIRST C COLUMN SPECIFIES THE VECTOR OF FUNCTION VALUES AND C THE SECOND THE VECTOR OF DERIVATIVES. C WORK - A WORKING STORAGE (DIMENSION IROW). C IROW - THE DIMENSION OF VECTORS Z AND WORK AND OF EACH C COLUMN IN MATRIX F. C ICOL - THE NUMBER OF COLUMNS IN F (I.E. 1 OR 2). C ARG - THE RESULTING VECTOR OF SELECTED AND ORDERED C ARGUMENT VALUES (DIMENSION NDIM). C VAL - THE RESULTING VECTOR OF SELECTED FUNCTION VALUES C (DIMENSION NDIM) IN CASE ICOL=1. IN CASE ICOL=2, C VAL IS THE VECTOR OF FUNCTION AND DERIVATIVE VALUES C (DIMENSION 2*NDIM) WHICH ARE STORED IN PAIRS (I.E. C EACH FUNCTION VALUE IS FOLLOWED BY ITS DERIVATIVE C VALUE). C NDIM - THE NUMBER OF POINTS WHICH MUST BE SELECTED OUT OF C THE GIVEN TABLE (Z,F). C C REMARKS C NO ACTION IN CASE IROW LESS THAN 1. C IF INPUT VALUE NDIM IS GREATER THAN IROW, THE PROGRAM C SELECTS ONLY A MAXIMUM TABLE OF IROW POINTS. THEREFORE THE C USER OUGHT TO CHECK CORRESPONDENCE BETWEEN TABLE (ARG,VAL) C AND ITS DIMENSION BY COMPARISON OF NDIM AND IROW, IN ORDER C TO GET CORRECT RESULTS IN FURTHER WORK WITH TABLE (ARG,VAL). C THIS TEST MAY BE DONE BEFORE OR AFTER CALLING C SUBROUTINE ATSG. C SUBROUTINE ATSG ESPECIALLY CAN BE USED FOR GENERATING THE C TABLE (ARG,VAL) NEEDED IN SUBROUTINES ALI, AHI, AND ACFI. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C SELECTION IS DONE BY GENERATING THE VECTOR WORK WITH C COMPONENTS WORK(I)=ABS(Z(I)-X) AND AT EACH OF THE NDIM STEPS C (OR IROW STEPS IF NDIM IS GREATER THAN IROW) C SEARCHING FOR THE SUBSCRIPT OF THE SMALLEST COMPONENT, WHICH C IS AFTERWARDS REPLACED BY A NUMBER GREATER THAN C MAX(WORK(I)). C C .................................................................. C SUBROUTINE ATSG(X,Z,F,WORK,IROW,ICOL,ARG,VAL,NDIM) C C DIMENSION Z(1),F(1),WORK(1),ARG(1),VAL(1) IF(IROW)11,11,1 1 N=NDIM C IF N IS GREATER THAN IROW, N IS SET EQUAL TO IROW. IF(N-IROW)3,3,2 2 N=IROW C C GENERATION OF VECTOR WORK AND COMPUTATION OF ITS GREATEST ELEMENT. 3 B=0. DO 5 I=1,IROW DELTA=ABS(Z(I)-X) IF(DELTA-B)5,5,4 4 B=DELTA 5 WORK(I)=DELTA C C GENERATION OF TABLE (ARG,VAL) B=B+1. DO 10 J=1,N DELTA=B DO 7 I=1,IROW IF(WORK(I)-DELTA)6,7,7 6 II=I DELTA=WORK(I) 7 CONTINUE ARG(J)=Z(II) IF(ICOL-1)8,9,8 8 VAL(2*J-1)=F(II) III=II+IROW VAL(2*J)=F(III) GOTO 10 9 VAL(J)=F(II) 10 WORK(II)=B 11 RETURN END C C .................................................................. C C SUBROUTINE ATSM C C PURPOSE C NDIM POINTS OF A GIVEN TABLE WITH MONOTONIC ARGUMENTS ARE C SELECTED AND ORDERED SUCH THAT C ABS(ARG(I)-X).GE.ABS(ARG(J)-X) IF I.GT.J. C C USAGE C CALL ATSM (X,Z,F,IROW,ICOL,ARG,VAL,NDIM) C C DESCRIPTION OF PARAMETERS C X - THE SEARCH ARGUMENT. C Z - THE VECTOR OF ARGUMENT VALUES (DIMENSION IROW). C THE ARGUMENT VALUES MUST BE STORED IN INCREASING C OR DECREASING SEQUENCE. C F - IN CASE ICOL=1, F IS THE VECTOR OF FUNCTION VALUES C (DIMENSION IROW). C IN CASE ICOL=2, F IS AN IROW BY 2 MATRIX. THE FIRST C COLUMN SPECIFIES THE VECTOR OF FUNCTION VALUES AND C THE SECOND THE VECTOR OF DERIVATIVES. C IROW - THE DIMENSION OF VECTOR Z AND OF EACH COLUMN C IN MATRIX F. C ICOL - THE NUMBER OF COLUMNS IN F (I.E. 1 OR 2). C ARG - THE RE