text:
      INCLUDE 'VICMAIN_FOR'
      SUBROUTINE MAIN44
C  IBIS ROUTINE IBISSTAT
C
C  PURPOSE:
C	IBISSTAT PERFORMS A VARIETY OF STATISTICAL ANALYSES
C		ON IBIS TABULAR FILES
C
C	ORIGINAL PROGRAMMER:    FRANK EVANS    OCTOBER 1986
C
C	VERSION 1:
C
C   1-87  SXP   Modified to use LLSQ to solve linear least
C               squares system of equations.
C   1-87  SXP   changed MAXTOTLEN to 200,000 to reduce memory usage
C               to prevent quota problems.
C   2-95  MAC   CRI  MSTP S/W CONVERSION (VICAR PORTING)
C
	IMPLICIT NONE
	INTEGER	MAXCOLLEN, MAXTOTLEN, MAXCOLS
	PARAMETER (MAXCOLLEN = 100000, MAXTOTLEN=200000, MAXCOLS=40)

	INTEGER	UNIT, OUTUNIT, STATUS, COUNT, DEF
	INTEGER	CLEN, NCOL, NINCOLS
	INTEGER	INCOLS(MAXCOLS), IN2COLS(MAXCOLS), COLS(MAXCOLS)
	INTEGER	I, II, J, JJ, K, M, N
	INTEGER	INDEX, BUFPTR
	INTEGER	NUMBINS, BIN, MAXCOUNTS, HISTCOUNT(1000)
	INTEGER	SCATSIZE(2), XWIDTH, YWIDTH
	INTEGER	XBIN, YBIN, SCATTER(75,50)
	INTEGER	SWATH, ROW, NUMCOL, STARTCOL
	INTEGER	DEPCOL, ERR, IPIV(MAXCOLS+1)
	INTEGER	ERRCOL, ERRDEF
	INTEGER	REPLICS(MAXCOLS), TOTREP, NDF(3)
	INTEGER	NA, NB, NUMPOINT(2)
	INTEGER	N1, N2, INIBIS, OUTIBIS, BEGIN

	REAL	DATA(MAXTOTLEN), DATA2(MAXTOTLEN)
	REAL	OUTDATA(MAXCOLLEN), RESID(MAXCOLLEN)
	REAL	DEPDATA(MAXCOLLEN), DEPDATA2(MAXCOLLEN)
	REAL	WORK(10000), DUMMY(10000)

	REAL	MEDIAN, MEAN, STANDEV, MINDATA, MAXDATA
	REAL	SCALE, OFFSET, BINVALUE(1000), HIST(1000), HISTFACTOR
	REAL	XSCALE, YSCALE, XOFFSET, YOFFSET
	REAL	MINX, MAXX, MINY, MAXY
	REAL	CORRARRAY(MAXCOLS*MAXCOLS), CORRMATRIX(MAXCOLS,MAXCOLS)
	REAL	CORRPROB(MAXCOLS,MAXCOLS), COVARRAY(MAXCOLS*MAXCOLS)
	REAL	R, DENOM, T
	REAL	SOLUTION(MAXCOLS+1), INTERVALS(MAXCOLS+1)
	REAL	SSR, RSQUARED, SEE
	REAL	DURWATSTAT, FRATIO, FPROB, PROB, CHISQR, CHIPROB
	REAL	EPSILON
	REAL	YBAR, EFFECTS(MAXCOLS), SSARRAY(3)
	REAL	T2, CRIT
	REAL	SUM
	REAL*8  EIGENCORR(MAXCOLS,MAXCOLS), EIGENVAL(MAXCOLS)
	REAL	PHI, CHI2, STUDIS, FISH

	LOGICAL NOPRINT, OUTPUT, XVPTST
	CHARACTER*8	OPTION, COLNAMES(MAXCOLS), DEPNAME, DIST
	CHARACTER*80	STRING, STARS, SCATLINE, SLINE

C		GET PARAMETERS
        CALL IFMESSAGE('IBISSTAT version 02-JAN-95')
	CALL XVP ('OUT', STRING, COUNT)
	OUTPUT = (COUNT .EQ. 1)

	NOPRINT = XVPTST ('NOPRINT')
	CALL XVP ('OPTION', OPTION, COUNT)
	CALL XVP ('COLS', INCOLS, NINCOLS)

	DO I = 1, MAXCOLS
	    IF (INCOLS(I) .LE. 9) THEN
		WRITE (COLNAMES(I), '(A6,I1)') 'COLUMN', INCOLS(I)
	    ELSE
		WRITE (COLNAMES(I), '(A6,I2)') 'COLUMN', INCOLS(I)
	    ENDIF
	ENDDO
	CALL XVPARM ('COLNAMES', COLNAMES, COUNT, DEF, 40)

	IF (.NOT. NOPRINT)   CALL XVMESSAGE (' ',' ')

C		OPEN THE INTERFACE FILE
        CALL XVUNIT(UNIT,'INP',1,STATUS,' ')
        IF (STATUS.NE.1) THEN
          CALL XVMESSAGE
     & (' INPUT FILE INITIAILIZATION ERROR-PROGRAM TERMINATED',' ')
          CALL ABEND
        END IF
        CALL IBIS_FILE_OPEN(UNIT,INIBIS,'READ',0,0,' ',' ',STATUS)
        IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(UNIT,STATUS,1)
        CALL IBIS_FILE_GET(INIBIS,'NR',CLEN,1,1)
        CALL IBIS_FILE_GET(INIBIS,'NC',NCOL,1,1)

	N = CLEN
	IF (N .LT. 2)  CALL MABEND (' COLUMN LENGTH TOO SHORT')
	IF (N .GT. MAXCOLLEN)  CALL MABEND (' COLUMN LENGTH TOO LONG')

	IF (OPTION(1:3) .EQ. 'SUM') GOTO 1000
	IF (OPTION(1:3) .EQ. 'HIS') GOTO 2000
	IF (OPTION(1:3) .EQ. 'SCA') GOTO 2500
	IF (OPTION(1:3) .EQ. 'COR') GOTO 3000
	IF (OPTION(1:3) .EQ. 'BEH') GOTO 4000
	IF (OPTION(1:3) .EQ. 'REG') GOTO 5000
	IF (OPTION(1:3) .EQ. 'ANO') GOTO 6000
	IF (OPTION(1:3) .EQ. 'FAC') GOTO 7000
	IF (OPTION(1:3) .EQ. 'DEN') GOTO 8000

C    *************************************************************************
C			STATISTICAL SUMMARY SECTION
1000	CONTINUE

	IF (OUTPUT) THEN
          CALL XVUNIT(OUTUNIT,'OUT',1,STATUS,' ')
          IF (STATUS.NE.1) THEN
            CALL XVMESSAGE
     &  (' OUTPUT FILE INITIAILIZATION ERROR-PROGRAM TERMINATED',' ')
            CALL ABEND
          END IF
          CALL IBIS_FILE_OPEN(OUTUNIT,OUTIBIS,'WRITE',NCOL,CLEN,' ',
     &           'COLUMN',STATUS)
          IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
	ENDIF

	IF (.NOT. NOPRINT) THEN
	    CALL XVMESSAGE ('              STATISTICAL SUMMARY ',' ')
	    WRITE (STRING, '(A,I6)') 'NUMBER OF CASES:', N
	    CALL XVMESSAGE (' ',' ')
	    CALL XVMESSAGE (STRING,' ')
            STRING(1:28) = 'COLUMN   NAME         MEDIAN'
            STRING(35:74) = 'MEAN      STD DEV        MIN         MAX'
	    CALL XVMESSAGE (STRING,' ')
	ENDIF

	DO INDEX = 1, NINCOLS
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 INCOLS(INDEX),STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,DATA,INCOLS(INDEX),
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)

	    CALL SORTR(DATA,N)

	    IF (MOD(N,2) .EQ. 0) THEN
		MEDIAN = ( DATA(N/2) + DATA(N/2+1) ) /2.0
	    ELSE
		MEDIAN = DATA((N+1)/2)
	    ENDIF
	    MINDATA = DATA(1)
	    MAXDATA = DATA(N)

	    MEAN = 0.0
	    DO I = 1, N
		MEAN = MEAN + DATA(I)
	    ENDDO
	    MEAN = MEAN/FLOAT(N)

	    STANDEV = 0.0
	    DO I = 1, N
		STANDEV = STANDEV + (DATA(I) - MEAN)**2
	    ENDDO
	    STANDEV = SQRT(STANDEV/(N-1.0))

	    IF (OUTPUT) THEN
		OUTDATA(1) = MEDIAN
		OUTDATA(2) = MEAN
		OUTDATA(3) = STANDEV
		OUTDATA(4) = MINDATA
		OUTDATA(5) = MAXDATA
                CALL IBIS_COLUMN_SET(OUTIBIS,'U_FORMAT',
     &                     'REAL',INDEX,STATUS)
                CALL IBIS_COLUMN_WRITE(OUTIBIS,OUTDATA,INDEX,
     &                     1,CLEN,STATUS)
                IF (STATUS.NE.1)
     &                     CALL IBIS_SIGNAL(OUTIBIS,STATUS,1)
	    ENDIF

	    IF (.NOT. NOPRINT) THEN
		WRITE (STRING, '(2X,I3,3X,A,5(1X,F11.4))')
     +			INCOLS(INDEX), COLNAMES(INDEX),
     +			MEDIAN, MEAN, STANDEV, MINDATA, MAXDATA
		CALL XVMESSAGE (STRING,' ')
	    ENDIF

	ENDDO

	IF (OUTPUT) THEN
          CALL IBIS_FILE_CLOSE(OUTIBIS,' ',STATUS)
          IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
        ENDIF

	GOTO 9999

C    *************************************************************************
C			HISTOGRAM SECTION
2000	CONTINUE

	CALL XVP ('BINS', NUMBINS, COUNT)

	IF (OUTPUT) THEN
          CALL XVUNIT(OUTUNIT,'OUT',1,STATUS,' ')
          IF (STATUS.NE.1) THEN
            CALL XVMESSAGE
     & (' OUTPUT FILE INITIAILIZATION ERROR-PROGRAM TERMINATED',' ')
            CALL ABEND
          END IF
          CALL IBIS_FILE_OPEN(OUTUNIT,OUTIBIS,'WRITE',NCOL,CLEN,' ',
     &           'COLUMN',STATUS)
          IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
        ENDIF

	IF (.NOT. NOPRINT) THEN
	    CALL XVMESSAGE ('      HISTOGRAM  ',' ')
	    CALL XVMESSAGE (' ',' ')
	    WRITE (STRING, '(A,I6)') 'NUMBER OF CASES:', N
	    CALL XVMESSAGE (STRING,' ')
	    CALL XVMESSAGE (' ',' ')
	ENDIF

	DO INDEX = 1, NINCOLS
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 INCOLS(INDEX),STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,DATA,INCOLS(INDEX),
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)

	    IF (.NOT. NOPRINT) THEN
		CALL XVMESSAGE (' ',' ')
		WRITE (STRING, '(A,A)')
     +			'       HISTOGRAM FOR : ', COLNAMES(INDEX)
		CALL XVMESSAGE (STRING,' ')
		CALL XVMESSAGE (' ',' ')
	    ENDIF

	    MINDATA = +1.0E38
	    MAXDATA = -1.0E38
	    DO I = 1, N
		MINDATA = MIN (MINDATA, DATA(I))
		MAXDATA = MAX (MAXDATA, DATA(I))
	    ENDDO
	    IF (MAXDATA .EQ. MINDATA) THEN
		WRITE (STRING, '(A,F10.3)')
     +			'       COLUMN MAX EQUALS MIN: ', MAXDATA
		CALL XVMESSAGE (STRING,' ')
		GOTO 2099
	    ENDIF
	    SCALE = (NUMBINS-0.01)/(MAXDATA-MINDATA)
	    OFFSET = -SCALE*MINDATA + 0.5

	    DO I = 1, N
		BIN = NINT(SCALE*DATA(I) + OFFSET)
		HISTCOUNT(BIN) = HISTCOUNT(BIN) + 1
	    ENDDO

	    DO BIN = 1, NUMBINS
		BINVALUE(BIN) =  (FLOAT(BIN) - OFFSET)/SCALE
		HIST(BIN) = FLOAT(HISTCOUNT(BIN))
		MAXCOUNTS = MAX (MAXCOUNTS,HISTCOUNT(BIN))
	    ENDDO

	    IF (OUTPUT) THEN
                CALL IBIS_COLUMN_SET(OUTIBIS,'U_FORMAT',
     &                     'REAL',2*INDEX-1,STATUS)
                CALL IBIS_COLUMN_WRITE(OUTIBIS,BINVALUE,2*INDEX-1,
     &                     1,NUMBINS,STATUS)
                IF (STATUS.NE.1)
     &                     CALL IBIS_SIGNAL(OUTIBIS,STATUS,1)
                CALL IBIS_COLUMN_SET(OUTIBIS,'U_FORMAT',
     &                     'REAL',2*INDEX,STATUS)
                CALL IBIS_COLUMN_WRITE(OUTIBIS,HIST,2*INDEX,
     &                     1,NUMBINS,STATUS)
                IF (STATUS.NE.1)
     &                     CALL IBIS_SIGNAL(OUTIBIS,STATUS,1)
	    ENDIF

	    IF (.NOT. NOPRINT) THEN
		STARS(1:1) = '|'
		HISTFACTOR = MIN (1.0, 60.0/MAXCOUNTS)
		DO I = 1, NUMBINS
	            STARS(2:79) = ' '
                    COUNT = NINT(HISTFACTOR*HISTCOUNT(I)) + 1
		    DO II = 2, COUNT
                      STARS(II:II) ='*'
                    ENDDO
                    WRITE (STRING, '(F10.3,1X,I5,2X,A60)')
     +				BINVALUE(I), HISTCOUNT(I), STARS
		    CALL XVMESSAGE (STRING,' ')
		ENDDO
		CALL XVMESSAGE (' ',' ')
	    ENDIF
2099	    CONTINUE
	ENDDO

	IF (OUTPUT) THEN
          CALL IBIS_FILE_CLOSE(OUTIBIS,' ',STATUS)
          IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
        ENDIF

	GOTO 9999

C    *************************************************************************
C			SCATTER PLOT SECTION
2500	CONTINUE

	IF (NOPRINT)  GOTO 9999

	CALL XVMESSAGE ('                  SCATTER PLOT  ',' ')
	CALL XVMESSAGE (' ',' ')
	WRITE (STRING, '(A,I6)') 'NUMBER OF CASES:', N
	CALL XVMESSAGE (STRING,' ')
	CALL XVMESSAGE (' ',' ')

	CALL XVP ('SCATSIZE', SCATSIZE, COUNT)
	XWIDTH = MIN (SCATSIZE(1), 75)
	YWIDTH = MIN (SCATSIZE(2), 50)

	DO INDEX = 1, NINCOLS-1, 2
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 INCOLS(INDEX),STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,DATA,INCOLS(INDEX),
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 INCOLS(INDEX+1),STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,DATA2,INCOLS(INDEX+1),
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)

	    WRITE (STRING, '(8X,A8,A,A8)')
     +			COLNAMES(INDEX), ' VS ', COLNAMES(INDEX+1)
	    CALL XVMESSAGE (STRING,' ')
	    CALL XVMESSAGE (' ',' ')

	    MINX = +1.0E38
	    MAXX = -1.0E38
	    DO I = 1, N
		MINX = MIN (MINX, DATA(I))
		MAXX = MAX (MAXX, DATA(I))
	    ENDDO
	    IF (MAXX .EQ. MINX) THEN
		CALL XVMESSAGE ('MAX X EQUALS MIN X',' ')
		GOTO 2599
	    ENDIF
	    XSCALE = (XWIDTH-0.01)/(MAXX-MINX)
	    XOFFSET = -XSCALE*MINX + 0.5

	    MINY = +1.0E38
	    MAXY = -1.0E38
	    DO I = 1, N
		MINY = MIN (MINY, DATA2(I))
		MAXY = MAX (MAXY, DATA2(I))
	    ENDDO
	    IF (MAXY .EQ. MINY) THEN
		CALL XVMESSAGE ('MAX Y EQUALS MIN Y',' ')
		GOTO 2599
	    ENDIF
	    YSCALE = (YWIDTH-0.01)/(MAXY-MINY)
	    YOFFSET = -YSCALE*MINY + 0.5

	    DO I = 1, N
		XBIN = NINT(XSCALE*DATA(I) + XOFFSET)
		YBIN = NINT(YSCALE*DATA2(I) + YOFFSET)
		SCATTER(XBIN,YBIN) = SCATTER(XBIN,YBIN) + 1
	    ENDDO

	    WRITE (STRING, '(F9.2)')  MAXY
	    CALL XVMESSAGE (STRING,' ')
	    DO YBIN = YWIDTH, 1, -1
		DO XBIN = 1, XWIDTH
		    COUNT = SCATTER(XBIN,YBIN)
		    COUNT = MIN(COUNT,9) + 48
		    IF (COUNT .EQ. 49)  COUNT = 42
		    IF (COUNT .EQ. 48)  COUNT = 32
		    SCATLINE(XBIN:XBIN) = CHAR(COUNT)
		ENDDO
		WRITE (STRING, '(A1,A60)') '+', SCATLINE
		CALL XVMESSAGE (STRING,' ')
	    ENDDO
	    STRING(1:1) = ' '
	    DO I = 1, XWIDTH+1
		STRING(I:I) = '+'
	    ENDDO
	    CALL XVMESSAGE(STRING,' ')
	    WRITE (STRING, '(F9.2)')  MINY
	    CALL XVMESSAGE(STRING,' ')
	    WRITE (STRING, '(8X,F9.2,30X,F9.2)')
     +			MINX, MAXX
	    CALL XVMESSAGE(STRING,' ')

2599	    CONTINUE
	    CALL XVMESSAGE (' ',' ')
	ENDDO

	GOTO 9999

C    *************************************************************************
C			CORRELATION SECTION
3000	CONTINUE

	IF (NINCOLS*N .GT. MAXTOTLEN)  CALL MABEND (' TOO MUCH DATA')

	DO J = 1, NINCOLS
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 INCOLS(J),STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,DATA((J-1)*N+1),INCOLS(J),
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)
	ENDDO
	M = NINCOLS

	CALL CORRE (N, M, 1, DATA,
     +			DUMMY(1), DUMMY(MAXCOLS+1), COVARRAY,
     +		  CORRARRAY, DUMMY(MAXCOLS**2+1), WORK, WORK(MAXCOLS+1) )

	K = 1
	DO I = 1, M
	    DO J = 1, I
		CORRMATRIX(I,J) = CORRARRAY(K)
		CORRMATRIX(J,I) = CORRMATRIX(I,J)
		K = K + 1
	    ENDDO
	ENDDO

	DO I = 1, M*M
	    COVARRAY(I) = COVARRAY(I)/N
	ENDDO

	DO I = 1, M
	    DO J = 1, I
		R = CORRMATRIX(I,J)
		DENOM = (1.0 - R**2) / (N-2)
		IF (DENOM .GT. 0.0) THEN
		    T = SQRT( R**2/DENOM )
		    CORRPROB(I,J) = 1.0 - STUDIS(N-2, T)
		ELSE
		    CORRPROB(I,J) = 0.0
		ENDIF
		CORRPROB(J,I) = CORRPROB(I,J)
	    ENDDO
	ENDDO

	IF (OUTPUT) THEN
	    DO I = 1, M
		COLS(I) = I
	    ENDDO
            CALL XVUNIT(OUTUNIT,'OUT',1,STATUS,' ')
            IF (STATUS.NE.1) THEN
            CALL XVMESSAGE
     & (' OUTPUT FILE INITIAILIZATION ERROR-PROGRAM TERMINATED',' ')
              CALL ABEND
            END IF
            CALL IBIS_FILE_OPEN(OUTUNIT,OUTIBIS,'WRITE',M,M,' ',
     &           'COLUMN',STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
	    DO I = 1, M
                CALL IBIS_COLUMN_SET(OUTIBIS,'U_FORMAT',
     &                     'REAL',I,STATUS)
                CALL IBIS_COLUMN_WRITE(OUTIBIS,COVARRAY(M*(I-1)+1),
     &                     I,1,M,STATUS)
                IF (STATUS.NE.1)
     &                     CALL IBIS_SIGNAL(OUTIBIS,STATUS,1)
	    ENDDO
            CALL IBIS_FILE_CLOSE(OUTIBIS,' ',STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
	ENDIF

	IF (.NOT. NOPRINT) THEN
	    CALL XVMESSAGE ('            CORRELATION   ',' ')
	    CALL XVMESSAGE (' ',' ')
	    WRITE (STRING, '(A,I6)') 'NUMBER OF CASES:', N
	    CALL XVMESSAGE (STRING,' ')
	    CALL XVMESSAGE (' ',' ')
	    DO SWATH = 1, INT((M-1)/6)+1
		STARTCOL = 6*(SWATH-1)+1
		NUMCOL = MIN(6, M-STARTCOL+1)
		CALL XVMESSAGE (' ',' ')

                STRING(1:80) = ' '
                SLINE(1:80) = ' '
                BEGIN = 1
                DO JJ = 1,NUMCOL
                  SLINE(BEGIN:BEGIN+7) = COLNAMES(JJ)
                  BEGIN = BEGIN + 11
                ENDDO

                STRING(14:10+(NUMCOL*11)) = SLINE
		CALL XVMESSAGE (STRING,' ')
		DO ROW = 1, M
		    K = MIN (ROW-STARTCOL+1,NUMCOL)
		    IF (K .GT. 0) THEN
                        STRING(1:80) = ' '
                        SLINE(1:80) = ' '
                        BEGIN = 12
                        STRING(1:8) = COLNAMES(ROW)
                        DO JJ = 1, K
                          WRITE (SLINE, '(2X,F7.4,2X)')
     +                           CORRMATRIX(JJ+STARTCOL-1,ROW)
                          STRING(BEGIN:BEGIN+11) = SLINE
                          BEGIN = BEGIN + 11
                        ENDDO
!			WRITE (STRING, '(3X,<K>(2X,F7.4,2X))')
!     +				(CORRMATRIX(J+STARTCOL-1,ROW), J = 1,K)
			CALL XVMESSAGE (STRING,' ')
                        STRING(1:80) = ' '
!			WRITE (STRING, '(11X,<K>(1X,F10.3))')
!     +			 (COVARRAY(J+STARTCOL-1+M*(ROW-1)), J = 1,K)
                        BEGIN = 12
                        DO JJ = 1, K
                          WRITE (SLINE, '(1X,F10.3)')
     +                           COVARRAY(JJ+STARTCOL-1+M*(ROW-1))
                          STRING(BEGIN:BEGIN+11) = SLINE
                          BEGIN = BEGIN + 11
                        ENDDO
			CALL XVMESSAGE (STRING,' ')
                        STRING(1:80) = ' '
                        BEGIN = 12
                        DO JJ = 1, K
                          WRITE (SLINE, '(2X,2HS=,F5.3,2X)')
     +                           CORRPROB(JJ+STARTCOL-1,ROW)
                          STRING(BEGIN:BEGIN+11) = SLINE
                          BEGIN = BEGIN + 11
                        ENDDO
!			WRITE (STRING, '(11X,<K>(2X,2HS=,F5.3,2X))')
!     +				(CORRPROB(J+STARTCOL-1,ROW), J = 1,K)
			CALL XVMESSAGE (STRING,' ')
			CALL XVMESSAGE (' ',' ')
		    ENDIF
		ENDDO
	    ENDDO
	ENDIF

	GOTO 9999

C    *************************************************************************
C			BEHRENS-FISHER TEST SECTION
4000	CONTINUE

	IF (NINCOLS*N .GT. MAXTOTLEN)  CALL MABEND (' TOO MUCH DATA')

	CALL XVP ('BCOLS', IN2COLS, COUNT)
	IF (COUNT .NE. NINCOLS)
     +		CALL MABEND (' TWO SAMPLES MUST HAVE SAME DIMENSION')

	CALL XVPARM ('NUMPOINT', NUMPOINT, COUNT, DEF, 2)
	IF (DEF .EQ. 0)  THEN
	    NA = NUMPOINT(1)
	    NB = NUMPOINT(2)
	ELSE
	    NA = N
	    NB = N
	ENDIF

	DO INDEX = 1, NINCOLS
	    BUFPTR = N*(INDEX-1) + 1
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 INCOLS(INDEX),STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,DATA(BUFPTR),INCOLS(INDEX),
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 IN2COLS(INDEX),STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,DATA2(BUFPTR),IN2COLS(INDEX),
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)
	ENDDO
	M = NINCOLS

	CALL BEHERN (DATA, NA, DATA2, NB, M, N,N, M,
     +		DUMMY, 0.01, T2, CORRARRAY,
     +		WORK(3*MAXCOLS+1), WORK(MAXCOLS+1), WORK(2*MAXCOLS+1),
     +		NDF, CRIT)

	FRATIO = T2*(NA-M)/(M*(NA-1))
	FPROB = FISH(FRATIO, M, NA-M)

	IF (.NOT. NOPRINT) THEN
	    CALL XVMESSAGE ('            BEHRENS-FISHER TEST  ',' ')
	    CALL XVMESSAGE (' ',' ')
	    WRITE (STRING, '(A,I6,1X,I6)') 'NUMBER OF DATA POINTS:',
     +						 NA, NB
	    CALL XVMESSAGE (STRING,' ')
	    CALL XVMESSAGE (' ',' ')
	    WRITE (STRING, '(A,F10.3,A,F7.3)')
     +			'HOTELLING T SQUARED = ', T2,
     +			'      PROB = ', FPROB
	    CALL XVMESSAGE (STRING,' ')
	ENDIF

	GOTO 9999

C    *************************************************************************
C			MULTIPLE REGRESSION SECTION
5000	CONTINUE

	IF ((NINCOLS+1)*N .GT. MAXTOTLEN)
     +			CALL MABEND (' TOO MUCH DATA TO FIT')
	IF (N .LT. NINCOLS+1)  CALL MABEND (' NOT ENOUGH ROWS OF DATA')

	DO I = 1, N
	    DATA(I) = 1.0
	ENDDO
	DO INDEX = 1, NINCOLS
	    BUFPTR = INDEX*N + 1
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 INCOLS(INDEX),STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,DATA(BUFPTR),INCOLS(INDEX),
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)
	ENDDO
	M = NINCOLS + 1

	CALL XVP ('DEPCOL', DEPCOL, COUNT)
	CALL XVP ('DEPNAME', DEPNAME, COUNT)
        CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 DEPCOL,STATUS)
        CALL IBIS_COLUMN_READ(INIBIS,DEPDATA,DEPCOL,
     &                 1,CLEN,STATUS)
        IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)

	DO I = 1, N*M
	    DATA2(I) = DATA(I)
	ENDDO
	DO I = 1, N
	    DEPDATA2(I) = DEPDATA(I)
	ENDDO

	EPSILON = 1.0E-7
	CALL LLSQ (DATA2, DEPDATA2, N, M, 1,
     +			SOLUTION, IPIV, EPSILON, ERR, WORK)
	IF (ERR .NE. 0)  CALL MABEND (' DEPENDENT COLUMNS')
	SSR = WORK(1)

	CALL CALCRES (DATA, DEPDATA, SOLUTION, RESID, N, M,
     +				SSR, DURWATSTAT)

	CALL REGSTAT (DATA, DEPDATA, N, M, SSR,
     +			SEE, RSQUARED, INTERVALS, FRATIO, FPROB)

	CALL XVPARM ('ERRCOL', ERRCOL, COUNT, ERRDEF, 1)
	IF (ERRDEF .EQ. 0) THEN
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 ERRCOL,STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,OUTDATA,ERRCOL,
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)
	    CHISQR = 0.0
	    DO I = 1, N
		CHISQR = CHISQR + ( RESID(I)/OUTDATA(I) )**2
	    ENDDO
	    CHIPROB = CHI2(CHISQR, FLOAT(N-M))
	    IF (CHIPROB .GT. 0.5)  CHIPROB = 1.0 - CHIPROB
	ENDIF

	PROB = 1.0 - (STUDIS(N-M, 1.0) - 0.5)

	IF (OUTPUT) THEN
	    DO I = 1, M
		OUTDATA(I) = SOLUTION(I)
	    ENDDO
          CALL XVUNIT(OUTUNIT,'OUT',1,STATUS,' ')
          IF (STATUS.NE.1) THEN
            CALL XVMESSAGE
     & (' OUTPUT FILE INITIAILIZATION ERROR-PROGRAM TERMINATED',' ')
            CALL ABEND
          END IF
          CALL IBIS_FILE_OPEN(OUTUNIT,OUTIBIS,'WRITE',NCOL,CLEN,' ',
     &           'COLUMN',STATUS)
          IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
                CALL IBIS_COLUMN_SET(OUTIBIS,'U_FORMAT',
     &                     'REAL',1,STATUS)
                CALL IBIS_COLUMN_WRITE(OUTIBIS,OUTDATA,1,1,
     &                     N,STATUS)
                IF (STATUS.NE.1)
     &                     CALL IBIS_SIGNAL(OUTIBIS,STATUS,1)
                CALL IBIS_COLUMN_SET(OUTIBIS,'U_FORMAT',
     &                     'REAL',2,STATUS)
                CALL IBIS_COLUMN_WRITE(OUTIBIS,RESID,2,
     &                     1,N,STATUS)
                IF (STATUS.NE.1)
     &                     CALL IBIS_SIGNAL(OUTIBIS,STATUS,1)
            CALL IBIS_FILE_CLOSE(OUTIBIS,' ',STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
	ENDIF

	IF (.NOT. NOPRINT) THEN
	    CALL XVMESSAGE ('            MULTIPLE REGRESSION  ',' ')
	    CALL XVMESSAGE (' ',' ')
	    WRITE (STRING, '(A,I6)') 'NUMBER OF CASES:', N
	    CALL XVMESSAGE (STRING,' ')
	    CALL XVMESSAGE (' ',' ')
	    WRITE (STRING, '(A,A)') 'DEPENDENT VARIABLE: ', DEPNAME
	    CALL XVMESSAGE (STRING,' ')
	    CALL XVMESSAGE (' ',' ')
	    WRITE (STRING, '(A,F7.3)')
     +			' VARIABLE     COEFFICIENT     ERROR: ', PROB
	    CALL XVMESSAGE (STRING,' ')
	    WRITE (STRING, '(2X,A,4X,F11.4,4X,F11.4)')
     +			'CONSTANT', SOLUTION(1), INTERVALS(1)
	    CALL XVMESSAGE (STRING,' ')
	    DO I = 2, M
		WRITE (STRING, '(2X,A,4X,F11.4,4X,F11.4)')
     +			COLNAMES(I-1), SOLUTION(I), INTERVALS(I)
		CALL XVMESSAGE (STRING,' ')
	    ENDDO
	    CALL XVMESSAGE (' ',' ')
	    WRITE (STRING, '(A,F7.4)') 'R SQUARED = ', RSQUARED
	    CALL XVMESSAGE (STRING,' ')
	    WRITE (STRING, '(A,F11.4)') 'STD ERR OF EST = ', SEE
	    CALL XVMESSAGE (STRING,' ')
	    WRITE (STRING, '(A,F11.4,A,F7.4)') 'OVERALL :   F RATIO = ',
     +			FRATIO, '      PROB = ', FPROB
	    CALL XVMESSAGE (STRING,' ')
	    WRITE (STRING, '(A,F7.4)')
     +			'DURBIN-WATSON STATISTIC = ', DURWATSTAT
	    CALL XVMESSAGE (STRING,' ')
	    IF (ERRDEF .EQ. 0) THEN
		WRITE (STRING, '(A,F12.4,A,F6.3)')
     +			'GOODNESS OF FIT:  CHI**2/DF = ', CHISQR/(N-M),
     +			'     PROB = ', CHIPROB
		CALL XVMESSAGE (STRING,' ')
	    ENDIF
	ENDIF

	GOTO 9999

C    *************************************************************************
C			ANOVA SECTION
6000	CONTINUE

	IF (NINCOLS*N .GT. MAXTOTLEN)  CALL MABEND (' TOO MUCH DATA')

	CALL XVPARM ('REPLICS', REPLICS, COUNT, DEF, 40)
	IF (DEF .EQ. 1) THEN
	    DO I = 1, NINCOLS
		REPLICS(I) = N
	    ENDDO
	ELSE IF (COUNT .NE. NINCOLS) THEN
	    CALL MABEND
     +		 (' MUST INPUT NUMBER OF REPLICATIONS FOR EACH GROUP')
	ENDIF

	BUFPTR = 1
	DO INDEX = 1, NINCOLS
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 INCOLS(INDEX),STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,DATA(BUFPTR),INCOLS(INDEX),
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)
	    BUFPTR = BUFPTR + REPLICS(INDEX)
	ENDDO
	TOTREP = BUFPTR - 1
	M = NINCOLS

	CALL ANOV1 (DATA, TOTREP, REPLICS, M, RESID, YBAR, EFFECTS,
     +			SSARRAY, NDF, DUMMY, FRATIO, FPROB, ERR)
	IF (ERR .NE. 0) CALL MABEND (' SUM OF SQUARES OVERFLOW')

	IF (.NOT. NOPRINT) THEN
	    CALL XVMESSAGE ('    ONEWAY ANALYSIS OF VARIANCE ',' ')
	    CALL XVMESSAGE (' ',' ')
	    WRITE (STRING, '(A,F11.4)') 'GRAND MEAN = ', YBAR
	    CALL XVMESSAGE (STRING,' ')
	    WRITE (STRING, '(A)')
     +			' VARIABLE     ESTIMATE OF EFFECT '
	    CALL XVMESSAGE (STRING,' ')
	    DO I = 1, M
		WRITE (STRING, '(2X,A,4X,F11.4)')
     +			COLNAMES(I), EFFECTS(I)
		CALL XVMESSAGE (STRING,' ')
	    ENDDO
	    CALL XVMESSAGE (' ',' ')
	    WRITE (STRING, '(A,F11.4,A,I4)')
     +		  'FACTOR:      SS = ', SSARRAY(1)
     +			,	'     DF = ', NDF(1)
	    CALL XVMESSAGE (STRING,' ')
	    WRITE (STRING, '(A,F11.4,A,I4)')
     +		  'RESIDUAL:    SS = ', SSARRAY(2)
     +			,	'     DF = ', NDF(2)
	    CALL XVMESSAGE (STRING,' ')
	    WRITE (STRING, '(A,F11.4,A,I4)')
     +		  'TOTAL:       SS = ', SSARRAY(3)
     +			,	'     DF = ', NDF(3)
	    CALL XVMESSAGE (STRING,' ')
	    WRITE (STRING, '(A,F11.4,A,F7.4)')
     +		'F RATIO = ', FRATIO
     +			, '      PROB = ', FPROB
	    CALL XVMESSAGE (STRING,' ')
	ENDIF

	GOTO 9999

C    *************************************************************************
C			FACTOR ANALYSIS SECTION
C			(PRINCIPAL COMPONENTS)
7000	CONTINUE

	IF (NINCOLS*N .GT. MAXTOTLEN)  CALL MABEND (' TOO MUCH DATA')

	DO J = 1, NINCOLS
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 INCOLS(J),STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,DATA((J-1)*N+1),INCOLS(J),
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)
	ENDDO
	M = NINCOLS

	CALL CORRE (N, M, 1, DATA,
     +			DUMMY(1), DUMMY(MAXCOLS+1), COVARRAY,
     +		  CORRARRAY, DUMMY(MAXCOLS**2+1), WORK, WORK(MAXCOLS+1) )

	K = 1
	DO I = 1, M
	    DO J = 1, M
		EIGENCORR(I,J) = DBLE(COVARRAY(K))/N
		K = K + 1
	    ENDDO
	ENDDO

	CALL SEIGEN (MAXCOLS, M, EIGENCORR, EIGENVAL,
     +					 WORK, EIGENCORR, ERR)
	IF (ERR .NE. 0) CALL MABEND (' DIAGONALIZATION FAILED')

	IF (.NOT. NOPRINT) THEN
	    CALL XVMESSAGE('   PRINCIPLE COMPONENTS FACTOR ANALYSIS  ',
     +				' ')
	    CALL XVMESSAGE (' ',' ')
	    CALL XVMESSAGE ('   FACTOR       EIGEN VALUE ',' ')
	    DO I = M, 1, -1
		WRITE (STRING, '(4X,I3,4X,F15.4)') M-I+1,SNGL(EIGENVAL(I))
		CALL XVMESSAGE (STRING,' ')
	    ENDDO
	    CALL XVMESSAGE (' ',' ')

	    DO SWATH = 1, INT((M-1)/6)+1
		STARTCOL = 6*(SWATH-1)+1
		NUMCOL = MIN(6, M-STARTCOL+1)
		CALL XVMESSAGE (' ',' ')
                STRING(1:80) = ' '
                SLINE(1:80) = ' '
                BEGIN = 1
                DO JJ = 1,NUMCOL
                  SLINE(BEGIN:BEGIN+7) = COLNAMES(JJ)
                  BEGIN = BEGIN + 11
                ENDDO
                STRING(15:10+(NUMCOL*11)) = SLINE
                CALL XVMESSAGE (STRING,' ')
                STRING(1:80) = ' '
                SLINE(1:80) = ' '

		DO J = 1, M
                  STRING(1:6) = 'FACTOR'
                  WRITE(SLINE, '(I2)') J
                  STRING(7:8) = SLINE
                  STRING(11:11) = '='
                  BEGIN = 15
                  DO JJ = 1, NUMCOL
                    WRITE(SLINE,'(F7.4)')
     +                SNGL(EIGENCORR(JJ+STARTCOL-1,M-J+1))
                    STRING(BEGIN:BEGIN+6) = SLINE
                    BEGIN = BEGIN + 11
                  ENDDO
!    MACOX      !!    WRITE (STRING, '(A6,I2,A4,<NUMCOL>(2X,F7.4,2X))')
!     +			'FACTOR', J, '  = ',
!     +			(SNGL(EIGENCORR(K+STARTCOL-1,M-J+1)), K = 1,NUMCOL)
		    CALL XVMESSAGE (STRING,' ')
		ENDDO
	    ENDDO
	ENDIF

	IF (OUTPUT) THEN
	    DO I = 1, N
		DO J = 1, M
		    SUM = 0.0
		    DO K = 1, M
			SUM = SUM + EIGENCORR(K,M-J+1)*DATA(N*(K-1)+I)
		    ENDDO
		    DATA2(N*(J-1)+I) = SUM
		ENDDO
	    ENDDO
          CALL XVUNIT(OUTUNIT,'OUT',1,STATUS,' ')
          IF (STATUS.NE.1) THEN
            CALL XVMESSAGE
     & (' OUTPUT FILE INITIAILIZATION ERROR-PROGRAM TERMINATED',' ')
            CALL ABEND
          END IF
          CALL IBIS_FILE_OPEN(OUTUNIT,OUTIBIS,'WRITE',NCOL,CLEN,' ',
     &           'COLUMN',STATUS)
          IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
	    DO J = 1, M
              CALL IBIS_COLUMN_SET(OUTIBIS,'U_FORMAT',
     &                     'REAL',J,STATUS)
              CALL IBIS_COLUMN_WRITE(OUTIBIS,DATA2(N*(J-1)+1),J,
     &                     1,N,STATUS)
              IF (STATUS.NE.1)
     &                     CALL IBIS_SIGNAL(OUTIBIS,STATUS,1)
	    ENDDO
            CALL IBIS_FILE_CLOSE(OUTIBIS,' ',STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
	ENDIF

	GOTO 9999

C    *************************************************************************
C			PROBABILITY DENSITY SECTION
8000	CONTINUE

	CALL XVP ('DISTRIB', DIST, COUNT)
	CALL XVP ('NDF', NDF, COUNT)

	IF (OUTPUT) THEN
          CALL XVUNIT(OUTUNIT,'OUT',1,STATUS,' ')
          IF (STATUS.NE.1) THEN
            CALL XVMESSAGE
     & (' OUTPUT FILE INITIAILIZATION ERROR-PROGRAM TERMINATED',' ')
            CALL ABEND
          END IF
          CALL IBIS_FILE_OPEN(OUTUNIT,OUTIBIS,'WRITE',2*NINCOLS,
     &           CLEN,' ','COLUMN',STATUS)
          IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
	ENDIF

	IF (.NOT. NOPRINT) THEN
	    CALL XVMESSAGE ('       PROBABILITY DENSITIES   ',' ')
	    CALL XVMESSAGE (' ',' ')
	ENDIF

	DO INDEX = 1, NINCOLS
            CALL IBIS_COLUMN_SET(INIBIS,'U_FORMAT','REAL',
     &                 INCOLS(INDEX),STATUS)
            CALL IBIS_COLUMN_READ(INIBIS,DATA,INCOLS(INDEX),
     &                 1,CLEN,STATUS)
            IF (STATUS.NE.1) CALL IBIS_SIGNAL(INIBIS,STATUS,1)

	    IF (DIST(1:3) .EQ. 'NOR') THEN
		DO I = 1, CLEN
		    DATA2(I) = PHI (DATA(I))
		ENDDO
		IF (.NOT. NOPRINT) THEN
		    WRITE (STRING, '(A)')
     +				'NORMAL (GAUSSIAN) DISTRIBUTION'
		ENDIF

	    ELSE IF (DIST(1:3) .EQ. 'CHI') THEN
		N1 = NDF(INDEX)
		DO I = 1, CLEN
		    DATA2(I) = CHI2 (DATA(I), FLOAT(N1) )
		ENDDO
		IF (.NOT. NOPRINT) THEN
		    WRITE (STRING, '(A,I3,A)')
     +				'CHI SQUARED DISTRIBUTION FOR ',
     +				N1, ' DEGREES OF FREEDOM'
		ENDIF

	    ELSE IF (DIST(1:3) .EQ. 'STU') THEN
		N1 = NDF(INDEX)
		DO I = 1, CLEN
		    DATA2(I) = STUDIS (N1, DATA(I) )
		ENDDO
		IF (.NOT. NOPRINT) THEN
		    WRITE (STRING, '(A,I3,A)')
     +				'STUDENTS T DISTRIBUTION FOR ',
     +				N1, ' DEGREES OF FREEDOM'
		ENDIF

	    ELSE IF (DIST(1:3) .EQ. 'FIS') THEN
		N1 = NDF(2*INDEX-1)
		N2 = NDF(2*INDEX)
		DO I = 1, CLEN
		    DATA2(I) = FISH (DATA(I), N1, N2)
		ENDDO
		IF (.NOT. NOPRINT) THEN
		    WRITE (STRING, '(A,I3,1X,I3,A)')
     +				'FISHERS F DISTRIBUTION FOR ',
     +				N1, N2, ' DEGREES OF FREEDOM'
		ENDIF
	    ENDIF

	    IF (OUTPUT) THEN
                CALL IBIS_COLUMN_SET(OUTIBIS,'U_FORMAT',
     &                     'REAL',2*INDEX-1,STATUS)
                CALL IBIS_COLUMN_WRITE(OUTIBIS,DATA,2*INDEX-1,
     &                     1,CLEN,STATUS)
                IF (STATUS.NE.1)
     &                     CALL IBIS_SIGNAL(OUTIBIS,STATUS,1)
                CALL IBIS_COLUMN_SET(OUTIBIS,'U_FORMAT',
     &                     'REAL',2*INDEX,STATUS)
                CALL IBIS_COLUMN_WRITE(OUTIBIS,DATA2,2*INDEX,
     &                     1,CLEN,STATUS)
                IF (STATUS.NE.1)
     &                     CALL IBIS_SIGNAL(OUTIBIS,STATUS,1)

	    ENDIF

	    IF (.NOT. NOPRINT) THEN
		CALL XVMESSAGE (STRING,' ')
	    ENDIF

	ENDDO

	IF (OUTPUT) THEN
          CALL IBIS_FILE_CLOSE(OUTIBIS,' ',STATUS)
          IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(OUTUNIT,STATUS,1)
        ENDIF

	GOTO 9999

9999	CONTINUE
        CALL IBIS_FILE_CLOSE(INIBIS,' ',STATUS)
        IF (STATUS.NE.1) CALL IBIS_SIGNAL_U(UNIT,STATUS,1)

	IF (.NOT. NOPRINT)   CALL XVMESSAGE(' ',' ')

	RETURN
	END

!	CHARACTER*16 FUNCTION FSTR (VALUE, FIELDWIDTH, NDECIMALS)
!	IMPLICIT NONE
!	REAL	VALUE
!	INTEGER	FIELDWIDTH, NDECIMALS
!
!	IF ( ABS(VALUE) .GE. 10**(FIELDWIDTH-NDECIMALS-2) .OR.
!     +	     (ABS(VALUE) .LT. 0.1**NDECIMALS .AND. VALUE .NE. 0.0) ) THEN
!	    WRITE (FSTR, '(E<FIELDWIDTH>.<NDECIMALS>)', ERR=100 )  VALUE
!	ELSE
!	    WRITE (FSTR, '(F<FIELDWIDTH>.<NDECIMALS>)', ERR=100 )  VALUE
!	ENDIF
!	RETURN
!
!100	FSTR = '****************'
!	RETURN
!	END

	SUBROUTINE REGSTAT(C, Y, N, M, SSR, SEE, RSQUARED, INTERVALS,
     +					FRATIO, FPROB)
	IMPLICIT NONE
	INTEGER	N, M
	REAL	C(N,M), Y(N), SSR, SEE, RSQUARED, INTERVALS(M)
	REAL	FRATIO, FPROB, FISH
	INTEGER	I, J, K, L, ERR, WORK2(41), WORK3(41)
	REAL	A(41*41), SUM, EPS, WORK1(41)
	REAL	MEAN, SUMSQR

	SUM = 0.0
	DO I = 1, N
	    SUM = SUM + Y(I)
	ENDDO
	MEAN = SUM/N
	SUMSQR = 0.0
	DO I = 1, N
	    SUMSQR = SUMSQR + (Y(I) - MEAN)**2
	ENDDO
	RSQUARED = (SUMSQR - SSR) / SUMSQR
	SEE = SQRT(SSR/(N-M))

C		CALCULATE CONFIDENCE INTERVALS FOR THE REGRESSION COEFS
	L = 1
	DO J = 1, M
	    DO K = 1, M
		SUM = 0.0
		DO I = 1, N
		    SUM = SUM + C(I,J)*C(I,K)
		ENDDO
		A(L) = SUM
		L = L + 1
	    ENDDO
	ENDDO

	EPS = 1.0E-7
	CALL MINV (A, M, EPS, ERR, WORK1, WORK2, WORK3)
	IF (ERR .NE. 0) CALL MABEND (' SINGULAR MATRIX')

	DO K = 1, M
	    INTERVALS(K) = SEE*SQRT(A(K+M*(K-1)))
	ENDDO

	IF (RSQUARED .EQ. 1.0) THEN
	    FRATIO = 99999999999.00
	    FPROB = 1.0
	ELSE
	    FRATIO = ( RSQUARED/(M-1) ) / ( (1-RSQUARED)/(N-M))
	    FPROB = FISH(FRATIO, M-1, N-M)
	ENDIF

	RETURN
	END

	SUBROUTINE CALCRES (A, B, X, RES, N, M, SSR, DURWATSTAT)
	IMPLICIT NONE
	INTEGER	I,J,  N,M
	REAL	A(N,M), B(N), X(M), RES(N), SSR, DURWATSTAT, SUM

	DO I = 1, N
	    SUM = B(I)
	    DO J = 1, M
		SUM = SUM - A(I,J)*X(J)
	    ENDDO
	    RES(I) = SUM
	ENDDO

	SSR = 0.0
	DO I = 1, N
	    SSR = SSR + RES(I)**2
	ENDDO
	SUM = 0.0
	DO I = 2, N
	    SUM = SUM + (RES(I)-RES(I-1))**2
	ENDDO
	DURWATSTAT = SUM/SSR

	RETURN
	END

      SUBROUTINE CORRE (N,M,IO,X,XBAR,STD,RX,R,B,D,T)
C
C ----PURPOSE:  COMPUTE MEANS, STANDARD DEVIATIONS, SUMS OF SQUARES
C     AND CROSS-PRODUCTS OF DEVIATIONS FROM MEANS, AND CORRELATION
C     COEFFICIENTS.
C
C ----CALLING SEQUENCE:  CALL CORRE(N,M,IO,X,XBAR,STD,RX,R,B,D,T)
C
C     N     NUMBER OF OBSERVATION.  N MUST BE GREATER THAN
C           OR EQUAL TO 2  (INPUT)
C
C     M     NUMBER OF VARIABLES  (INPUT)
C
C     IO    =0, IF A USER PROVIDED SUBROUTINE NAMED DATA(M,D) IS TO BE
C           USED TO READ THE DATA FROM A SPECIFIC INPUT DEVICE
C           =1, IF ALL DATA ARE IN CORRE IN X  (INPUT)
C
C     X     THE N BY M DATA MATRIX (INPUT IF IO=1)
C
C     XBAR  THE VECTOR OF LENGTH M CONTAINING MEANS  (OUTPUT)
C
C     STD   THE VECTOR OF LENGTH M CONTAINING STANDARD DEVIATIONS
C           (OUTPUT)
C
C     RX    THE M BY M MATRIX CONTAINING SUMS OF SQUARES AND CROSS-
C           PRODUCTS OF DEVIATIONS FROM MEANS. (SSCP MATRIX) (OUTPUT)
C
C     R     OUTPUT VECTOR CONTAINING CORRELATION COEFFICIENTS.
C           ONLY THE UPPER TRIANGULAR PORTION OF THE CORRELATION
C           MATRIX IS STORED.  THE COEFFICIENT IN THE I-TH ROW AND
C           THE J-TH COLUMN IS IN POSITION K=(I*(I-1)/2)+J OF R(K).
C
C     B     THE VECTOR OF LENGTH M CONTAINING THE DIAGONAL ELEMENTS
C           OF R, I.E. THE SUMS OF SQUARES   (OUTPUT)
C
C     D,T   WORKING VECTORS OF LENGTH M
C
C ----RESTRICTIONS AND COMMENTS:  THE USER MUST SUPPLY A
C     SUBROUTINE CALLED DATA(M,D)
C
C   --IF IO=0, DATA(M,D) READS ONE VECTOR, D, OF LENGTH M,
C     FROM A SPECIFIC INPUT DEVICE, I.E. ONE OF N OBSERVATION VECTORS,
C     EACH OF LENGTH M, IS STORED IN D EACH TIME DATA(M,D) IS CALLED.
C
C   --IF IO=1, AN UNUSED SUBROUTINE NAMED DATA MUST BE SUPPLIED.
C     FOR EXAMPLE:
C           SUBROUTINE DATA
C           RETURN
C           END
      DIMENSION X(1),XBAR(1),STD(1),RX(1),R(1),B(1),D(1),T(1)
C
C ----INITIALIZATION
      DO 100 J=1,M
      B(J)=0.0
  100 T(J)=0.0
      K=(M*M+M)/2
      DO 102 I=1,K
  102 R(I)=0.0
      FN=N
      L=0
C
C ----TEST IO
      IF(IO.EQ.0) GO TO 127
C
C ----DATA ARE ALREADY IN CORE
      DO 108 J=1,M
      DO 107 I=1,N
      L=L+1
      T(J)=T(J)+X(L)
  107 CONTINUE
      XBAR(J)=T(J)
      T(J)=T(J)/FN
  108 CONTINUE
C
      DO 115 I=1,N
      JK=0
      L=I-N
      DO 110 J=1,M
      L=L+N
      D(J)=X(L)-T(J)
  110 B(J)=B(J)+D(J)
      DO 115 J=1,M
      DO 115 K=1,J
      JK=JK+1
      R(JK)=R(JK)+D(J)*D(K)
  115 CONTINUE
      GO TO 205
C
C ----STORE PART OF THE DATA IN RX AND COMPUTE PARTIAL MEANS
  127 IF(N-M) 130, 130, 135
  130 KK=N
      GO TO 137
  135 KK=M
  137 DO 140 I=1,KK
C      CALL DATA (M,D)
      DO 140 J=1,M
      T(J)=T(J)+D(J)
      L=L+1
  140 RX(L)=D(J)
      FKK=KK
      DO 150 J=1,M
      XBAR(J)=T(J)
  150 T(J)=T(J)/FKK
C
C ----CALCULATE SUMS OF CROSS-PRODUCTS OF DEVIATIONS
C     FROM TEMPORARY MEANS FOR M OBSERVATIONS
      L=0
      DO 180 I=1,KK
      JK=0
      DO 170 J=1,M
      L=L+1
  170 D(J)=RX(L)-T(J)
      DO 180 J=1,M
      B(J)=B(J)+D(J)
      DO 180 K=1,J
      JK=JK+1
      R(JK)=R(JK)+D(J)*D(K)
  180 CONTINUE
C
C ----TEST KK
      IF(N-KK) 205, 205, 185
C
C ----READ THE REST OF OBSERVATIONS ONE AT A TIME, SUM
C     THE OBSERVATION, AND CALCULATE SUMS OF CROSS-
C     RPODUCTS OF DEVIATIONS FROM TEMPORARY MEANS
  185 KK=N-KK
      DO 200 I=1,KK
      JK=0
C      CALL DATA (M,D)
      DO 190 J=1,M
      XBAR(J)=XBAR(J)+D(J)
      D(J)=D(J)-T(J)
  190 B(J)=B(J)+D(J)
      DO 200 J=1,M
      DO 200 K=1,J
      JK=JK+1
      R(JK)=R(JK)+D(J)*D(K)
  200 CONTINUE
C
C ----CALCULATE MEANS
  205 JK=0
      DO 210 J=1,M
      XBAR(J)=XBAR(J)/FN
C
C ----ADJUST SUMS OF CROSS-PRODUCTS OF DEVIATIONS
C     FROM TEMPORARY MEANS
      DO 210 K=1,J
      JK=JK+1
      R(JK)=R(JK)-B(J)*B(K)/FN
  210 CONTINUE
C
C ----CALCULATE CORRELATION COEFFICIENTS
      JK=0
      DO 220 J=1,M
      JK=JK+J
  220 STD(J)= SQRT( ABS(R(JK)))
      DO 230 J=1,M
      DO 230 K=J,M
      JK=J+(K*K-K)/2
      L=M*(J-1)+K
      RX(L)=R(JK)
      L=M*(K-1)+J
      RX(L)=R(JK)
      IF(STD(J)*STD(K)) 225, 222, 225
  222 R(JK)=0.0
      GO TO 230
  225 CONTINUE
      R(JK)=R(JK)/(STD(J)*STD(K))
  230 CONTINUE
C
C ----CALCULATE STANDARD DEVIATIONS
      FN=SQRT(FN-1.0)
      DO 240 J=1,M
  240 STD(J)=STD(J)/FN
C
C ----COPY THE DIAGONAL OF THE MATRIX OF SUMS OF CROSS-PRODUCTS OF
C     DEVIATIONS FROM MEANS
      L=-M
      DO 250 I=1,M
      L=L+M+1
  250 B(I)=RX(L)
      RETURN
      END

      SUBROUTINE ANOV1(Y,N,J,I,R,YBAR,S,SS,NDF,SM,F,PF,IND)
C
C ----PURPOSE:  PERFORM ONE WAY ANALYSIS OF VARIANCE WITH AN UNEQUAL
C     NUMBER OF REPLICATIONS IN EACH CELL
C
C ----CALLING SEQUENCE: ANOV1(Y,N,J,I,R,YBAR,S,SS,NDF,SM,F,PF,IND)
C
C     Y     IS THE ARRAY OF DATA
C
C     N     IS THE TOTAL NUMBER OF REPLICATIONS
C
C     J     IS THE ARRAY OF THE NUMBER OF REPLICATIONS FOR EACH GROUP
C
C     I     IS THE NUMBER OF GROUPS
C
C     R     IS THE ARRAY OF RESIDUALS
C
C     YBAR  IS THE GRAND MEAN
C
C     S     IS THE ESTIMATE OF EFFECTS ARRAY
C
C     SS    IS THE SUMS OF SQUARES ARRAY
C
C     NDF   IS THE ARRAY OF THE DEGREES OF FREEDOM FOR EACH SUM
C           OF SQUARES
C
C     SM    IS THE  MEAN SQUARES ARRAY
C
C     F     IS THE F-STATISTIC
C
C     PF    IS THE PROBABILITY THAT AN F-DISTRIBUTED VARIABLE EXCEEDS F
C
C     IND   =0, IF OVERFLOW DOES NOT OCCUR
C           =1, IF OVERFLOW COMPUTING THE SUMS OF SQUARES OCCURS
C
C ----RESTRICTIONS AND COMMENTS: FISH AND ERFF ARE REQUIRED
      DIMENSION Y(N),R(N),J(I),SS(3),NDF(3),SM(3) ,S(I)
      XMAX=1.E30
      IND=1
C
C ----INITIALIZATION
      SS(1)=0.0
      SS(2)=0.0
      SS(3)=0.0
      LJ=1
      MJ=0
C
C ----CALCULATE THE GRAND MEAN
      YBAR=0.0
      DO 3 KI=1,N
    3 YBAR=YBAR+Y(KI)
      YBAR=YBAR/N
C
C ----CALCULATE THE GROUP MEANS
      DO 2 KI=1,I
      MJ=MJ+J(KI)
      YX=0.0
      DO 4  K=LJ,MJ
    4 YX=YX+Y( K)
      YX=YX/J(KI)
C
C ----CALCULATE THE GROUP EFFECTS
      S(KI)=YX-YBAR
C
C ----CALCULATE THE GROUP SUM OF SQUARES
      SS(1)=SS(1)+S(KI)*S(KI)*J(KI)
      DO 1 K=LJ,MJ
      R(K)=Y(K)-YX
C
C ----CALCULATE RESIDUAL SUM OF SQUARES
      SS(2)=SS(2)+R(K)*R(K)
C
C ----CALCULATE TOTAL SUM OF SQUARES
      SS(3)=SS(3)+(Y(K)-YBAR)*(Y(K)-YBAR)
      IF(SS(3).GT.XMAX) RETURN
    1 CONTINUE
    2 LJ=MJ+1
C
C ----CALCULATE THE DEGREES OF FREEDOM
      NDF(1)=I-1
      NDF(2)=N-I
      NDF(3)=N-1
C
C ----CALCULATE THE MEAN SQUARES
      SM(1)=SS(1)/NDF(1)
      SM(2)=SS(2)/NDF(2)
C
C ----CALCULATE THE F-STATISTIC
      F=SM(1)/SM(2)
C
C ----CALCULATE THE PROBABILITY THAT THE F-STATISTIC IS EXCEEDED
C     UNDER THE HYPOTHESIS THAT THERE IS NO DIFFERENCE BETWEEN THE
C     GROUPS
      PF=FISH(F,NDF(1),NDF(2))
      IND=0
      RETURN
      END

      SUBROUTINE BEHERN(A,NA,B,NB,M,NNA,NNB,MM,X,AL,T2,S,U,AF,AT,N2,
     .                  CRIT)
C
C ----PURPOSE: BEHERN  IS A SUBROUTINE FOR TESTING THE EQUALITY OF
C     TWO MULTIVARIATE NORMAL MEANS WHEN COVARIANCE MATRICES ARE
C     ASSUMED UNEQUAL.
C
C ----CALLING SEQUENCE: CALL BEHERN(A,NA,B,NB,M,NNA,NNB,MM,X,AL,T2,
C     S,U,AF,AT,N2,CRIT)
C
C     A     THE NA BY M ARRAY OF DATA FROM THE FIRST POPULATION
C
C     NA    THE NUMBER OF DATA POINTS IN THE FIRST POPULATION
C
C     B     THE NB BY M ARRAY OF DATA FROM THE SECOND POPULATION
C
C     NB    THE NUMBER OF DATA POINTS IN THE SECOND POPULATION
C
C     M     THE NUMBER OF VARIABLES (DIMENSIONS) IN THE
C           MULTIVARIATE POPULATIONS
C
C     NNA   THE MAXIMUM FIRST DIMENSION OF A
C
C     NNB   THE MAXIMUM FIRST DIMENSION OF B
C
C     MM   THE MAXIMUM FIRST DIMENSION OF S AND U
C
C     X     THE VECTOR OF DIMENSION M WHICH ARE THE VALUES OF THE
C           NULL HYPOTHESIS OF MEAN DIFFERENCES
C
C     AL    THE CONFIDENCE INTERVAL TYPE 1 ERROR
C
C     T2    THE HOTELLING T**2 VALUE
C
C     S     THE M BY M WORK ARRAY CONTAINING THE INVERSE OF THE
C           VARIANCE-COVARIANCE MATRIX, ON RETURN
C
C     U     AN M BY M WORK ARRAY
C
C     AF    IS AN M DIMENSIONAL WORK VECTOR
C
C     AT    AN M DIMENSIONAL WORK VECTOR
C
C     N2    DEGREES OF FREEDOM FOR TESTING (OUTPUT)
C
C     CRIT  THE COMPUTED CRITICAL VALUE FOR THE GIVEN AL (OUTPUT)
C
C ----RESTRICTIONS AND COMMENTS:  ONE MAY NOT HAVE PRIOR KNOWLEDGE THAT
C     COVARIANCE-A.NE.COVARIANCE-B  (SEE SUBROUTINE T2DIFF FOR EQUAL
C     COVARIANCE ASSUMPTION)
      REAL*8 S,U,AF,AT
      DIMENSION A(NNA,MM),B(NNB,MM),X(MM),S(MM,MM),U(MM,MM),AF(MM),
     2          AT(MM)
C
C ----CALCULATE VARIANCE-COVARIANCE MATRICES
      CALL COVARI(A,NA,M,NNA,MM,S,AF)
      CALL COVARI(B,NB,M,NNB,MM,U,AT)
C
C ----POOL VARIANCE-COVARIANCE MATRICES
      DO 30 I=1,M
   30 X(I)=X(I)+AF(I)-AT(I)
      MIN=NA
      IF(MIN.LE.NB) GO TO 21
      MIN=NB
      CALL POOLT2(B,NB,A,NA,M,NNB,NNA,MM,X,S,AT,AF,MIN,U)
      GO TO 29
   21 CALL POOLT2(A,NA,B,NB,M,NNA,NNB,MM,X,S,AF,AT,MIN,U)
   29 CONTINUE
C
      TOL=-1.
      CALL GINV (S,MM,MM,M,M,K,U,AF,AT)
C
      DO 40 I=1,M
      AF(I)=0.
      DO 40 J=1,M
   40 AF(I)=AF(I) + X(J)*S(J,I)
      T2=0.
      DO 50 I=1,M
   50 T2=T2 + AF(I)*X(I)
      T2=MIN*T2
      N2=MIN-M
      IF(N2.GT.0) GO TO 58
      RETURN
   58 CRIT= FLOAT((MIN-1)*M)*FISHIN(AL,M,N2)/N2
      RETURN
      END

      SUBROUTINE COVARI(A,N,M,NN,MM,S,X)
C
C ----PURPOSE: TO COMPUTE A COVARIANCE AS REQUIRED BY BEHER AND T2DIFF
C
C ----CALLING SEQUENCE: CALL COVARI(A,N,M,NN,MM,S,X)
C
C     A     THE NN BY MM ARRAY HOLDING THE N BY M ARRAY
C           OF N DATA VECTORS EACH OF LENGTH M  (INPUT)
C
C     N     THE NUMBER OF OBSERVED VECTORS  (INPUT)
C
C     M     THE LENGTH OF EACH OBSERVATION VECTOR  (INPUT)
C
C     NN    THE MAXIMUM NUMBER OF OBSERVATIONS  (INPUT)
C     MM    THE MAXIMUM LENGTH OF THE OBSERVATION VECTORS (INPUT)
C     S     THE MM BY MM ARRAY HOLDING THE M BY M COVARIANCE MATRIX
C           (OUTPUT)
C
C     X     THE MM VECTOR HOLDING THE M SAMPLE MEANS  (OUTPUT)
C
C ----RESTRICTIONS AND COMMENTS:  NONE
      REAL*8 S,X
      DIMENSION A(NN,MM),S(MM,MM),X(MM)
C
C ----INITIALIZATION
      MN=M-1
      DO 1 I=1,MN
      L=I+1
      DO 1 J=L,M
    1 S(I,J)=0.0
      DO 3 I=1,M
      SS=0.0
C
C ----COMPUTE SAMPLE MEANS
      DO 2 K=1,N
    2 SS=SS+A(K,I)
    3 X(I)=SS/FLOAT(N)
C
C ----COMPUTE THE COVARIANCE MATRIX
      DO 5 I=1,M
      SS=0.0
      DO 4 K=1,N
    4 SS=SS+(A(K,I)-X(I))**2
    5 S(I,I)=SS
      DO 7 K=1,N
      L=0
      DO 6 I=1,MN
      L=I+1
      DO 6 J=L,M
    6 S(I,J)=S(I,J)+(A(K,I)-X(I))*(A(K,J)-X(J))
    7 CONTINUE
      DO 8 I=2,M
      L=I-1
      DO 8 J=1,L
    8 S(I,J)=S(J,I)
      DO 9 I=1,M
      DO 9 J=1,M
    9 S(I,J)=S(I,J)/FLOAT(N-1)
      RETURN
      END

      SUBROUTINE POOLT2(A,NA,B,NB,M,NNA,NNB,MM,X,S,AF,AT,MIN,WORK)
C
C ----PURPOSE: TO COMPUTE THE POOLED COVARINACE MATRIX FOR THE BEHERN
C     TEST WHEN THE COVARIANCE MATRICES ARE ASSUMED UNEQUAL.
C
C ----CALLING SEQUENCE: CALL POOLT2(A,NA,B,NB,M,NNA,NNB,MM,X,S,AF,AT,MIN
C                                    ,WORK)
C
C     A     THE NA BY M ARRAY OF DATA FOR THE FIRST POPULATION (INPUT)
C
C     NA    THE NUMBER OF DATA POINTS OF THE FIRST POPULATION (NA LESS
C           THAN OR EQUAL TO NB) (INPUT)
C
C     B     THE NB BY M ARRAY OF DATA FOR THE SECOND POPULATION (INPUT)
C
C     NB    THE NUMBER OF DATA POINTS OF THE SECOND POPULATION (NB
C           GREATER THAN OR EQUAL TO NA) (INPUT)
C
C     M     THE NUMBER OF VARIABLES (DIMENSION) IN THE MULTIVARIATE
C           POPULATION (INPUT)
C
C     NNA   THE MAXIMUM FIRST DIMENSION OF A (INPUT)
C
C     NNB   THE MAXIMUM FIRST DIMENSION OF B (INPUT)
C
C     MM    THE MAXIMUM FIRST DIMENSION OF S AND U (INPUT)
C
C     X     THE WORKING ARRAY OF MEANS: X(I)=X(I)+AF(I)-AT(I), WHERE THE
C           X(I) ON THE RIGHT ARE CONSTANTS SPECIFIED BY THE NULL
C           HYPOTHESIS. (INPUT)
C
C     S     THE M BY M POOLED COVARIANCE MATRIX RETURNED TO THE USER
C           (OUTPUT)
C
C     AF    THE VECTOR OF SAMPLE MEANS FOR THE FIRST POPULATION (INPUT)
C
C     AT    THE VECTOR OF SAMPLE MEANS FOR THE SECOND POPULATION (INPUT)
C
C     MIN   THE MINIMUM OF NA AND NB (INPUT)
C
C     WORK  A WORKING ARRAY OF DIMENSION M
C
C ----RESTRICTIONS AND COMMENTS:THE ARRAYS S, AF, AT, AND WORK ARE
C     DOUBLE PRECISION ARRAYS.
      REAL*8 S,AF,AT,WORK
      DIMENSION A(NNA,MM),B(NNB,MM),X(MM),S(MM,MM),AF(MM),AT(MM),
     *WORK(MM)
C
      SNBNA=1./(MIN-1.)
      IF(NA.EQ.NB) GO TO 100
      SNA=1./NA
      SNAOB=SQRT(FLOAT(NA)/FLOAT(NB))
C
      DO 10 I=1,M
      W=0.
      DO 9 J=1,MIN
    9 W=W+B(J,I)
   10 WORK(I)= W*SNA
      DO 20 I=1,M
      DO 20 J=I,M
      W=0.
      DO 19 K=1,MIN
   19 W=W+(A(K,I)-AF(I)-X(I)-SNAOB*(B(K,I)-WORK(I)))*
     2    (A(K,J)-AF(J)-X(J)-SNAOB*(B(K,J)-WORK(J)))
      S(I,J)=W*SNBNA
   20 S(J,I)=S(I,J)
      RETURN
  100 DO 120 I=1,M
      DO 120 J=I,M
      W=0.
      DO 119 K=1,NA
  119 W=W+(A(K,I)-AF(I)-X(I)-B(K,I)+AT(I))*
     2    (A(K,J)-AF(J)-X(J)-B(K,J)+AT(J))
      S(I,J)=W*SNBNA
  120 S(J,I)=S(I,J)
      RETURN
      END

       FUNCTION STUDIS(N,T)
C
C ----PURPOSE: TO COMPUTE AN APPROXIMATION TO THE CUMMULATIVE DISTRIBU-
C     TION FUNCTION OF STUDENT'S T-DISTRIBUTION
C
C ----CALLING SEQUENCE: VAR=STUDIS(N,T)
C
C     N     THE DEGREES OF FREEDOM,N.GE.1  (INPUT)
C
C     T     THE VALUE AT WHICH THE CUMULATIVE PROBABILITY IS DESIRED
C           (INPUT)
C
C    STUDIS THE VALUE OF THE FUNCTION (OUTPUT)
C
C ----RESTRICTIONS AND COMMENTS;  FOR LARGE N USE NORMAL APPROX-
C     MATIONS
C
C ----INITIALIZATION
C
       AN=N
       TH=ATAN(T/SQRT(AN))
C
C ----TESTN = 1
C
       W=COS(TH)
       A=W*W
      IF(N.EQ.1) GOTO 60
      J=N-2
C
C ----TEST N EVEN
C
      IF (MOD(N,2).EQ.0) GOTO13
C
C ----COMPUTE EXPAN FOR N ODD AND N.GE.3
C
       SU=W
      IF (N.EQ.3) GOTO 11
       R=0.
       V=1.
       DO 10 I=3,J,2
       R=R+2.0
       V=V+2.0
       W=W*(R/V)*A
       SU=SU+W
   10 CONTINUE
   11 A=2./3.1415927*(TH+SIN(TH)*SU)
       GOTO 50
C
C ----COMPUTE STUDIS FOR EVEN
C
   13  W=1.
       SU=1.
      IF (N.EQ.2) GOTO 40
       R=-1.
       V=0.
       DO 20 I=2,J,2
       R=R+2.
       V=V+2.
       W=W*R/V*A
       SU=W+SU
   20  CONTINUE
   40  A=SIN(TH)*SU
       GOTO 50
   60 A=2./3.1415927*TH
C
C ----CHANGE FROM PR(ABS(X).LE.T) TO PR(X.LE.T)
C
   50  CONTINUE
      A=0.5*A+0.5
      STUDIS=A
       RETURN
       END

      FUNCTION CHI2(X,V)
C
C-----PURPOSE:  TO COMPUTE AN APPROXIMATE VALUE OF THE CHI-SQUARE
C     CUMULATIVE DISTRIBUTION FUNCTION GIVEN A POINT AND THE
C     DEGREES OF FREEDOM.
C
C-----CALLING SEQUENCE:  VARIABLE = CHI2(X,V)
C
C     X     THE POINT AT WHICH THE FUNCTION IS TO BE EVALUATED.
C     V     DEGREES OF FREEDOM  (REAL,POSITIVE)
C     CHI2  COMPUTED VALUE OF THE FUNCTION.
C ----RESTRICTIONS AND COMMENTS:  FUNCTIONS GAMINC AND PHI ARE REQUIRED
C
      IF (X.LT.160.) GO TO 5
      Z = (X-V)/SQRT(2.*V)
      CHI2 = PHI(Z)
      RETURN
C
    5 X2=X/2.0
      A=V/2.0
      CHI2=GAMINC(X2,A)
      RETURN
      END

      FUNCTION GAMINC(X,A)
C
C ----PURPOSE:  TO COMPUTE AN APPROXIMATION TO THE CUMULATIVE
C     DISTRIBUTION FUNCTION OF A GAMMA DISTRIBUTION, I.E.
C     APPROXIMATE INCOMPLETE GAMMA INTEGRALS.
C
C ----CALLING SEQUENCE:  VARIABLE = GAMINC(X,A)
C
C     X     THE VALUE OF THE RANDOM VARIABLE, X.GE.0  (INPUT)
C
C     A     THE VALUE OF THE PARAMETER, A.GE.0  (INPUT)
C
C     GAMINC THE VALUE OF THE CUMULATIVE DISTRIBUTION FUNCTION (OUTPUT)
C
C ----RESTRICTIONS AND COMMENTS:  NONE
      DATA A1,A2,A3,A4,A5,A6,A7,A8/-.5771916  ,.9882058  ,-.8970569  ,
     1.9182068  ,-.7567040  ,.4821993  ,-.1935278  ,.0358683  /
      DATA RELEPS/1.E-7/
      IF(X.GT.0.) GO TO 1
      GAMINC=0.0
      RETURN
    1 N=A
      F=A-N
C
C ----COMPUTE (X**F)EXP(-X)/GAMMA(F+1) FOR COEFFICIENTS
C     SEE EQN. 6.1.36,P. 257  ABRAMOWITZ AND STEGUN
      P=EXP(-X+F*ALOG(X))/(1.+F*(A1+F*(A2+F*(A3+F*(A4+F*(A5+F*(A6+F*(A7+
     1F*A8))))))))
C
C ----COMPUTE P*(1+F)...(N+F) TO COMPLETE THE COMPUTATION OF GAMMA (A+1)
      IF(N.EQ.0)  GO TO 5
      DO 2 I=1,N
    2 P=P*X/(F+I)
C
C ----INITIALIZE GAMINC
    5 GAMINC=P
      IF (P.EQ.0.) GO TO 6
      F=A+1.
C
C ----COMPUTE NEXT TERM IN THE APPROXIMATING SUM
    3 P=P*X/F
C
C ----TEST RELATIVE ERROR
      IF((P/GAMINC).LT.RELEPS) GO TO 4
C
C ----ACCUMULATE SUM
      GAMINC=GAMINC+P
      F=F+1.0
      GO TO 3
    4 RETURN
    6 GAMINC = -1.0
      RETURN
      END

      FUNCTION FISH(F,NN1,NN2)
C
C ----PURPOSE: TO COMPUTE THE CUMULATIVE DISTRIBUTION FUNCTION FOR
C     FISHER F-DISTRIBUTION
C
C ----CALLING SEQUENCE:  VARIABLE = FISH(F,N1,N2)
C     F     THE RATION (U/N1)/(V/N2), WHERE U AND V ARE INDEPENDENT
C           CHI-SQUARE RANDOM VARIABLES WITH N1 AND N2 DEGREES
C           OF FREEDOM, RESPECTIVELY.  (INPUT)
C
C     N1    THE FIRST DEGREES OF FREEDOM, N1.GE.1  (INPUT)
C
C     N2    THE SECOND DEGREES OF FREEDOM, N2.GE.1  (INPUT)
C
C ----RESTRICTIONS AND COMMENTS:  SUBROUTINE ERFF IS REQUIRED TO
C     EVALUATE PHI(X)
      LOGICAL E1,E2,E3
      N1 = NN1
      N2 = NN2
       IF(N1.GE.100.AND.N2.GE.100) GOTO 9
C
C ----INITIALIZATION AND SETTING OF LOGICAL SWITCHES TO .TRUE. IF
C     THE DEGREES OF FREEDOM ARE EVEN
      E1=.FALSE.
      E2=.FALSE.
      E3=.FALSE.
      IF(MOD(N1,2).EQ.0) E1=.TRUE.
      IF(MOD(N2,2).EQ.0) E2=.TRUE.
      X=N2/(N2+N1*F)
      IF(.NOT.(E1.OR.E2)) GO TO 5
      IF(E1.AND..NOT.E2) GO TO 1
      IF(.NOT.E1.AND.E2) GO TO 2
      IF(N1.LE.N2) GO TO 1
C
C ----INITIALIZATION FOR SECOND DEGREE OF FREEDOM EVEN AND LESS THAN
C     FIRST DEGREE OF FREEDOM IF IT TOO IS EVEN
    2 I=N1
      N1=N2
      N2=I
      X=1.0-X
      E3=.TRUE.
C
C ----INITIALIZATION FOR FIRST DEGREE OF FREEDOM EVEN AND LESS THAN
C     SECOND DEGREE OF FREEDOM IF IT IS EVEN
    1 Y=1.0-X
C
C ----CALCULATION OF PROBABILITY FOR AT LEAST ONE DEGREE OF FREEDOM
C     EVEN
      FISH=0.0
      H=SQRT(X**N2)
      M=N1/2
      DO  3 I=1,M
      FISH=FISH+H
    3 H=(H*Y*(N2+2.*(I-1))) / (2.*I)
      IF(E3) GO TO 4
C
C ----ADJUST CALCULATED PROBABILITY IF ITS ONES COMPLEMENT WAS
C     CALCULATED ORIGINALLY
      FISH=1.0-FISH
      RETURN
    4 I=N1
      N1=N2
      N2=I
      RETURN
C
C ----CALCULATION OF THE PROBABILITY FOR BOTH DEGREES OF FREEDOM ODD
    5 Y=1.0-X
      H=.63661977*SQRT(X*Y)
      FISH=.63661977* ACOS(SQRT(X))
      IF(N2.EQ.1) GO TO 8
      M=N2-2
      DO  6  I=1,M,2
      FISH=FISH+H
    6 H=H*X*(I+1)/(I+2)
    8 IF(N1.EQ.1) RETURN
      H=H*N2
      M=N1-2
      DO  7  I=1,M,2
      FISH=FISH-H
    7 H=H*Y*(N2+I)/(I+2)
      RETURN
    9  D1=N1
       D2=N2
       DT=(D1/D2)*F
       DN=SQRT((2.*D2-1.)*DT)-SQRT(2.*D1-1.)
       X=DN/SQRT(1.+DT)
       FISH=PHI(X)
       RETURN
      END

      FUNCTION ERFF(X,N1)
C ----PURPOSE:  TO COMPUTE THE VALUE OF THE ERROR FUNCTION OR
C     THE STANDARD NORMAL CUMULATIVE DISTRIBUTION FUNCTION
C     AT A GIVEN POINT.
C
C ----CALLING SEQUENCE :  VARIABLE = ERFF(X,N1)
C
C     X     IS THE POINT WHERE THE FUNCTION IS EVALUATED (INPUT)
C
C     N1    =1, IF THE ERROR FUNCTION IS DESIRED.
C           =2, IF THE COMPLIMENTARY ERROR FUNCTION IS DESIRED.
C           =3, IF THE STANDARD NORMAL CUMULATIVE DISTRIBUTION
C           FUNCTION IS DESIRED.
C
C ----RESTRICTIONS AND COMMENTS: NONE.
      REAL ABS,T,V
      GO TO (1,2,3),N1
C
C ----ERF
      ENTRY ERF(X)
    1 N=1
      T=X
      GO TO 4
C
C ----ERFC
      ENTRY ERFC(X)
    2 N=2
      T=X
      GO TO 4
C
C ----PHI
      ENTRY PHI(X)
    3 N=3
      T=X/1.4142136
    4 M=1
      IF(X.LT.0.) M=2
      V=ABS(T)
      U=T*T
C
C ----TEST V
      IF(V-0.5) 10,10,20
C
C ----CASE V.LE.0.5
   10 ERFF= ((0.31665289*U+1.7222758)*U+21.385332)/((U+7.8437457)*U
     1    +18.952257)*T
C
C ----BRANCH TO RETURN THE APPROPRIATE FUNCTION
      GO TO (5,6,7),N
    5 RETURN
    6 ERFF=1.0-ERFF
      RETURN
    7 ERFF=0.5+0.5*ERFF
      RETURN
C
   20 X1=V
      X2=V*V
      X3=V*X2
      X4=X2*X2
C
C ----TEST V
      IF(V-8.0) 25,25,30
C
C ----CASE  V BETWEEN 0.5 AND 8.0 (HART #5703,P 294)
   25 ERFF=EXP(-X2)*(6.1337529+X1*6.1772458+X2*2.8501393+X3*0.56409092)
     1  /  (6.1337546+X1*13.098327+X2*11.497651+X3*5.0472398+X4)
      GO TO 40
C
C ----EXCEPT FOR THE LAST CASE, ERFF=0
   30 ERFF=0.
C
C ----TEST N AND V
      IF((N.NE.2.AND.M.EQ.1).OR.(N.NE.3.AND.M.EQ.2).OR.(V.GT.13.))GOTO40
C
C ----CASE  V BETWEEN 8. AND 13.5 AND ERFC OR PHI(-V) IS DESIRED
C     (HART #5721,P 296)
      ERFF=EXP(-X2)*(0.14845921+X1*0.56418774)/(0.51143728+
     1     X1*0.26277059+X2)
C
C ----BRANCH IF X.LT.0 ,OTHERWISE BRANCH TO RETURN THE
C     APPROPRIATE FUNCTION.
   40 IF(M.GT.1) GO TO 50
      GO TO(41,5,42),N
   41 ERFF=1.0-ERFF
      RETURN
   42 ERFF=0.5-0.5*ERFF+0.5
      RETURN
C
C ----BRANCH TO RETURN THE APPROPRIATE FUNCTION
   50 GO TO (51,52,53),N
   51 ERFF=ERFF-1.0
      ERF=ERFF
      RETURN
   52 ERFF=2.0-ERFF
      ERFC=ERFF
      RETURN
   53 ERFF=0.5*ERFF
      PHI=ERFF
      RETURN
      END

       FUNCTION FISHIN(ALPHA,N1,N2)
C
C ----PURPOSE:  TO COMPUTE THE INVERSE OF FISHER'S
C     F-DISTRIBUTION, GIVEN THE CONFIDENCE COEFFICIENT.
C
C ----CALLING SEQUENCE:  VARIABLE = FISHIN (ALPHA,N1,N2)
C
C     FISHIN  THE NAME OF THE FUNCTION CONTAINING THE
C             VALUE OF THE INVERSE (OUTPUT)
C
C     ALPHA   THE SIGNIFICANCE LEVEL (E.G. .01,.05.,ETC) (INPUT)
C
C     N1      THE DEGREES OF FREEDOM OF THE NUMERATOR (INPUT)
C
C     N2      THE DEGREES OF FREEDOM OF THE DENOMINATOR (INPUT)
C ----RESTRICTIONS AND COMMENTS:  THE SUBROUTINES FISH, PHINV, AND
C     ERFF ARE REQUIRED.
C
C ----CALCULATES THE INVERSE 'F' VALUE GIVEN THE CONFIDENCE COEFFICIENT
C     ALPHA AND THE DEGREES OF FREEDOM(N).
      Y1=N1
      Y2=N2
C
C ----ADJUST FOR DEGREES OF FREEDOM EQUAL TO 1
      IF(N1.EQ.1)  Y1=2
      IF(N2.EQ.1)  Y2=2
C
C ----CALL PHINV TO GET INVERSE NORMAL VALUE OF 1.-ALHPA
       X=PHINV(1.-ALPHA)
C
C ----COMPUTE LAMDA VALUE
      Y=(X**2-3.)/6.
      IC=0
C
C ----COMPUTE THE INITIAL APPROXIMATION TO THE INVERSE 'F' FUNCTION
      Y1=1./(Y1-1.)
      Y2=1./(Y2-1.)
      H=2./(Y1+Y2)
      X=X*SQRT(H+Y)/H-(Y1-Y2)*(Y+5./6.-2./(3.*H))
      X=EXP(2.*X)
      ISIZ=1
      IF(N1.GT.100.OR.N2.GT.100) ISIZ=2
C
C ----COMPUTE THE CONSTANT TO THE 'F' DISTRIBUTION, TESTING FOR N1 AND/
C     OR N2 ODD OR EVEN.
      G=1.
      IB1=2
      IF(MOD(N1,2).EQ.0) GO TO 1
      G=1.7724539
      IB1=1
    1 IB2=2
      IF(MOD(N2,2).EQ.0) GO TO 2
      G=G*1.7724539
      IB2=1
    2 IB3=2
      IF(MOD(N1+N2,2).EQ.0) GO TO 3
      G=G/1.7724539
      IB3=1
    3 IF((IB1+IB2).NE.2)  G=2.*G
      IF((N1+N2).LE.3)  GO TO 5
      ND=N1+N2-2-IB3 + 1
      DO 4 J=1,ND,2
      I=J-1
      IF((IB1+I).LE.(N1-2))  G=G*(IB1+I)
      IF((IB2+I).LE.(N2-2))  G=G*(IB2+I)
    4 G=G/(IB3+I)
C
C ----COMPUTE THE VALUE OF FISHIN
      GO TO (5,55),ISIZ
    5 Y2=N2/(N2+N1*X)
      Y1=1.-Y2
      Y=1.+(G*(1.-ALPHA-FISH(X,N1,N2)))/SQRT(Y1**N1*Y2**N2)
      GO TO 6
 55   Y2=N2/(N2+N1*X)
      Y1=(N1*X)/(N2+N1*X)
      Z=((X**(1./3.)*((9.*N2-2.)/(9.*N2)))-((9.*N1-2.)/(9.*N1)))/SQRT((2
     *./9.)*(1./N1+(X**(2./3.))*(1./N2)))
      F=PHI(Z)
      DENOM=N1/2.*ALOG(Y1)+N2/2.*ALOG(Y2)
      Y=1.+EXP(ALOG(G)-DENOM)*(1.0-ALPHA-F)
    6 FISHIN=X*Y
C
C ----IF FISHIN IS NEGATIVE, RESET FISHIN TO .5*LAST APPROXIMATION(X).
      IF(Y.LT.0.)  FISHIN=.5*X
C
C ----IF THE ABSOLUTE VALUE OF THE DIFFERENCE IS LESS THAN .5E-6, RETURN
      IF(ABS(X/FISHIN-1.).LT.(.5E-6))  GO TO 7
C
C ----IF THE RELATIVE VALUE OF THE DIFFERENCE IS LESS THAN .5E-6, RETURN
      IF(ABS(X-FISHIN).LT.(.5E-6)) GO TO 7
      IC=IC+1
      IF(IC.GT.100)  RETURN
C
C ----SET THE APPROXIMATION EQUAL TO FISHNN AND CONTINUE TO ITERATE
      X=FISHIN
      GO TO (5,55),ISIZ
    7 RETURN
      END

      FUNCTION PHINV( PR )
C
C ----PURPOSE:  TO COMPUTE THE INVERSE OF THE STANDARD NORMAL CUMULATIVE
C     DISTRIBUTION FUNCTION (MEAN 0, VARIANCE 1)
C
C ----CALLING SEQUENCE:  VARIABLE=PHINV(P)
C     P     THE PROBABILITY FOR WHICH THE INVERSE
C           IS TO BE CALCULATED (INPUT)
C
C -----RESTRICTIONS AND COMMENTS:  THE SUBPROGRAM ERFF(PHI)
C     IS REQUIRED
C
C ----TEST P
      P = PR
      IF (P.GT.0.0. AND.P.LT.1.0) GO TO 5
      IF (P.EQ.1.0) PHINV=1.0E+38
      IF(P.EQ.0.0) PHINV=-1.0E+38
C
C ----ERROR IF P NOT BETWEEN ZERO AND ONE
      IF(P.LT.0.0.OR.P.GT.1.0) PHINV=1.0E20
      RETURN
C
C ----SET FLAG:  X POSITIVE
C     WHEN K=1, NEGATIVE
C     WHEN K=2.
    5 K = 1
      IF(P .GT. 0.5)GO TO 47
C
C ----COMPUTE FIRST APPROXIMATION, EQUATION 26.2.23, PAGE 933,
C     ABRAMOWITZ AND STEGUN
    8 T3=SQRT(-2.0*ALOG(P))
      T4P=2.515517+.802853*T3+.010328*T3*T3
      T5P=1.0+1.432788*T3+.189269*T3*T3+.001308*T3*T3*T3
      XT=T3-T4P/T5P
      XT=-XT
C
C ----NEWTON RHAPSON
C     ITERATION
   13 DO 53 I=1,100
      PHP = EXP(-0.5*XT*XT)
      PT = PHI (XT)
      IF(ABS(P-PT) .LT. P*4.0E-8)GO TO 99
      Z = (P-PT)*2.50662827 / PHP
      XT = XT + Z
   53 CONTINUE
      GO TO 99
   47  P = 1.0 - P
      K = 2
      GO TO 8
C
C ----ASSIGN SIGN
C     + OR - TO X.
   99 GO TO (26,27),K
   26 PHINV = XT
      RETURN
   27 PHINV = -XT
      P = 1.0 - P
      RETURN
      END

      SUBROUTINE SEIGEN(NM,N,A,D,E,Z,IERR)
C
      INTEGER I,J,K,L,M,N,II,JP1,NM,MML
      REAL*8 A(NM,N),D(N),E(N),Z(NM,N)
      REAL*8 B,C,F,G,H,P,R,S,HH,SCALE,MACHEP
      REAL*8 DSQRT,DABS,DSIGN
C
C     THIS SUBROUTINE IS A MERGER OF THE EISPACK SUBROUTINES TRED2 AND
C     TQL2, WHICH ARE THEMSELVES TRANSLATED FROM THE ALGOL PROCEDURES.
C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C     MATRIX EIGENSYSTEM ROUTINES - EISPACK GUIDE, 30(1974).
C
C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS OF A
C     FULL SYMMETRIC REAL MATRIX.
C
C     THE WORK IS DONE IN TWO PARTS.  THE FIRST (CORRESPONDING TO
C     TRED2) REDUCES THE FULL SYMMETRIC MATRIX TO A SYMMETRIC
C     TRIDIAGONAL MATRIX USING AND ACCUMULATING ORTHOGONAL SIMILARITY
C     TRANSFORMATIONS.  THE SECOND (TQL2) FINDS THE EIGENVALUES AND
C     EIGENVECTORS OF THE SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
C
C     ON INPUT:
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT;
C
C        N IS THE ORDER OF THE MATRIX;
C
C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C        E IS A REAL VECTOR FOR INTERMEDIATE STORAGE.
C
C     ON OUTPUT:
C
C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
C          UNORDERED FOR INDICES 1,2,...,IERR-1;
C
C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
C          EIGENVALUES;
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C
C        A AND Z MAY COINCIDE.  IF DISTINCT, A IS UNALTERED.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     ------------------------------------------------------------------
C
      DO 100 I = 1, N
C
         DO 100 J = 1, I
            Z(I,J) = A(I,J)
  100 CONTINUE
C
      IF (N .EQ. 1) GO TO 320
C     :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- ::::::::::
      DO 300 II = 2, N
         I = N + 2 - II
         L = I - 1
         H = 0.0D0
         SCALE = 0.0D0
         IF (L .LT. 2) GO TO 130
C     :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) ::::::::::
         DO 120 K = 1, L
  120    SCALE = SCALE + DABS(Z(I,K))
C
         IF (SCALE .NE. 0.0D0) GO TO 140
  130    E(I) = Z(I,L)
         GO TO 290
C
  140    DO 150 K = 1, L
            Z(I,K) = Z(I,K) / SCALE
            H = H + Z(I,K) * Z(I,K)
  150    CONTINUE
C
         F = Z(I,L)
         G = -DSIGN(DSQRT(H),F)
         E(I) = SCALE * G
         H = H - F * G
         Z(I,L) = F - G
         F = 0.0D0
C
         DO 240 J = 1, L
            Z(J,I) = Z(I,J) / (SCALE * H)
            G = 0.0D0
C     :::::::::: FORM ELEMENT OF A*U ::::::::::
            DO 180 K = 1, J
  180       G = G + Z(J,K) * Z(I,K)
C
            JP1 = J + 1
            IF (L .LT. JP1) GO TO 220
C
            DO 200 K = JP1, L
  200       G = G + Z(K,J) * Z(I,K)
C     :::::::::: FORM ELEMENT OF P ::::::::::
  220       E(J) = G / H
            F = F + E(J) * Z(I,J)
  240    CONTINUE
C
         HH = F / (H + H)
C     :::::::::: FORM REDUCED A ::::::::::
         DO 260 J = 1, L
            F = Z(I,J)
            G = E(J) - HH * F
            E(J) = G
C
            DO 260 K = 1, J
               Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K)
  260    CONTINUE
C
         DO 280 K = 1, L
  280    Z(I,K) = SCALE * Z(I,K)
C
  290    D(I) = H
  300 CONTINUE
C
  320 D(1) = 0.0D0
      E(1) = 0.0D0
C     :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES ::::::::::
      DO 500 I = 1, N
         L = I - 1
         IF (D(I) .EQ. 0.0D0) GO TO 380
C
         DO 360 J = 1, L
            G = 0.0D0
C
            DO 340 K = 1, L
  340       G = G + Z(I,K) * Z(K,J)
C
            DO 360 K = 1, L
               Z(K,J) = Z(K,J) - G * Z(K,I)
  360    CONTINUE
C
  380    D(I) = Z(I,I)
         Z(I,I) = 1.0D0
         IF (L .LT. 1) GO TO 500
C
         DO 400 J = 1, L
            Z(I,J) = 0.0D0
            Z(J,I) = 0.0D0
  400    CONTINUE
C
  500 CONTINUE
C
C     :::::::::: LAST CARD OF TRED2 ::::::::::
C
C     ------------------------------------------------------------------
C
C     :::::::::: MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C                IT IS THE SMALLEST POSITIVE NUMBER SUCH THAT 1.+MACHEP=
C                MACHEP = 16.0D0**(-13) FOR LONG FORM ARITHMETIC
C                ON S360 ::::::::::
C                MACHEP = 2.0D0**(-55) FOR LONG FORM ARITHMETIC ON PDP
      MACHEP=2.775557564D-17
C
      IERR = 0
      IF (N .EQ. 1) GO TO 1001
C
      DO1100 I = 2, N
 1100 E(I-1) = E(I)
C
      F = 0.0D0
      B = 0.0D0
      E(N) = 0.0D0
C
      DO1240 L = 1, N
         J = 0
         H = MACHEP * (DABS(D(L)) + DABS(E(L)))
         IF (B .LT. H) B = H
C     :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT ::::::::::
         DO1110 M = L, N
            IF (DABS(E(M)) .LE. B) GO TO1120
C     :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP ::::::::::
 1110    CONTINUE
C
 1120    IF (M .EQ. L) GO TO1220
 1130    IF (J .EQ. 30) GO TO 1000
         J = J + 1
C     :::::::::: FORM SHIFT ::::::::::
         P = (D(L+1) - D(L)) / (2.0D0 * E(L))
         R = DSQRT(P*P+1.0D0)
         H = D(L) - E(L) / (P + DSIGN(R,P))
C
         DO1140 I = L, N
 1140    D(I) = D(I) - H
C
         F = F + H
C     :::::::::: QL TRANSFORMATION ::::::::::
         P = D(M)
         C = 1.0D0
         S = 0.0D0
         MML = M - L
C     :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- ::::::::::
         DO1200 II = 1, MML
            I = M - II
            G = C * E(I)
            H = C * P
            IF (DABS(P) .LT. DABS(E(I))) GO TO1150
            C = E(I) / P
            R = DSQRT(C*C+1.0D0)
            E(I+1) = S * P * R
            S = C / R
            C = 1.0D0 / R
            GO TO1160
 1150       C = P / E(I)
            R = DSQRT(C*C+1.0D0)
            E(I+1) = S * E(I) * R
            S = 1.0D0 / R
            C = C * S
 1160       P = C * D(I) - S * G
            D(I+1) = H + S * (C * G + S * D(I))
C     :::::::::: FORM VECTOR ::::::::::
            DO1180 K = 1, N
               H = Z(K,I+1)
               Z(K,I+1) = S * Z(K,I) + C * H
               Z(K,I) = C * Z(K,I) - S * H
 1180       CONTINUE
C
 1200    CONTINUE
C
         E(L) = S * P
         D(L) = C * P
         IF (DABS(E(L)) .GT. B) GO TO1130
 1220    D(L) = D(L) + F
 1240 CONTINUE
C     :::::::::: ORDER EIGENVALUES AND EIGENVECTORS ::::::::::
      DO1300 II = 2, N
         I = II - 1
         K = I
         P = D(I)
C
         DO1260 J = II, N
            IF (D(J) .GE. P) GO TO1260
            K = J
            P = D(J)
 1260    CONTINUE
C
         IF (K .EQ. I) GO TO1300
         D(K) = D(I)
         D(I) = P
C
         DO1280 J = 1, N
            P = Z(J,I)
            Z(J,I) = Z(J,K)
            Z(J,K) = P
 1280    CONTINUE
C
 1300 CONTINUE
C
      GO TO 1001
C     :::::::::: SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS ::::::::::
 1000 IERR = L
 1001 RETURN
C     :::::::::: LAST CARD OF TQL2 ::::::::::
C
C     :::::::::: LAST CARD OF SEIGEN ::::::::::
      END

      SUBROUTINE MINV  (A,N,E,K,X,J2,I2)
C*     MATRIX INVERSION ROUTINE-FORMULATED BY E. G. CLAYTON
C
C      A--SQUARE ARRAY (SINGLE PRECISION) CONTAINING ORIGINAL MATRIX
C      N--ORDER OF ORIGINAL MATRIX
C      E--TEST CRITERION FOR NEAR ZERO DIVISOR
C      K--LOCATION FOR SINGULARITY OR ILL-CONDITION INDICATOR
C         K=0 =) MATRIX NONSINGULAR.
C         K=1 =) MATRIX SINGULAR (OR ILL-CONDITIONED)
C      X--A WORK VECTOR OF SIZE N
C      J2--AN INTEGER WORK VECTOR OF SIZE N
C      I2--AN INTEGER WORK VECTOR OF SIZE N
C
      DIMENSION A(N,N),X(N),J2(N),I2(N)
C
C     INITIALIZATION
C
      K=1
      I2(1)=0
      J2(1)=0
C
C     BEGIN COMPUTATION OF INVERSE
C
      DO 15 L=1,N
      L1=L-1
      IF(L1.EQ.0)L1=1
      BIGA=-1.0
C
C     LOOK FOR THE ELEMENT OF GREATEST ABSOLUTE VALUE,CHOOSING
C     ONE FROM A ROW AND COLUMN NOT PREVIOUSLY USED.
C
      DO 5 I=1,N
      DO 1 I3=1,L1
      IF(I .EQ. I2(I3))  GO TO 5
    1 CONTINUE
      DO 4 J=1,N
      DO 2 I3=1,L1
      IF(J .EQ. J2(I3))  GO TO 4
    2 CONTINUE
      AT=ABS(A(I,J))
      IF(BIGA .GE. AT)  GO TO 4
      BIGA=AT
      J1=J
      I1=I
    4 CONTINUE
    5 CONTINUE
C
C     TAG THE ROW AND COLUMN FROM WHICH THE ELEMENT IS CHOSEN.
C
      J2(L)=J1
      I2(L)=I1
      DIV=A(I1,J1)
C
C     TEST ELEMENT AGAINST ZERO CRITERION
      IF(ABS(DIV) .LE. E)  GO TO 221
C
C     PERFORM THE COMPUTATIONS
C
      OOD = 1./DIV
      DO 7 J=1,N
      A(I1,J)=A(I1,J)*OOD
    7 CONTINUE
      A(I1,J1)=OOD
      DO 11 I=1,N
      IF(I1 .EQ. I)  GO TO 11
      DO 10 J=1,N
      IF(J1 .EQ. J)  GO TO 10
      A(I,J)=A(I,J)-A(I1,J)*A(I,J1)
   10 CONTINUE
   11 CONTINUE
      DO 14 I=1,N
      IF(I1 .EQ. I)  GO TO 14
      A(I,J1)=-A(I,J1)*A(I1,J1)
   14 CONTINUE
   15 CONTINUE
C
C     COMPUTATION COMPLETE AT THIS POINT
C     UNSCRAMBLE THE INVERSE
C
      DO 18 J=1,N
      DO 16 I=1,N
      I1=I2(I)
      J1=J2(I)
      X(J1)=A(I1,J)
   16 CONTINUE
      DO 17 I=1,N
      A(I,J)=X(I)
17    CONTINUE
   18 CONTINUE
      DO 21 I=1,N
      DO 19 J=1,N
      I1=I2(J)
      J1=J2(J)
      X(I1)=A(I,J1)
   19 CONTINUE
      DO 20 J=1,N
      A(I,J)=X(J)
   20 CONTINUE
   21 CONTINUE
      K=0
  221 RETURN
      END

      SUBROUTINE GINV  (A,MR,MC,NR,NC,KZ,U,AFLAG,ATEMP)
C-----------------------------------------------------------------------
C
C          THIS ROUTINE COMPUTES THE GENERALIZED INVERSE OF A MATRIX
C          AND STORES THE TRANSPOSE OF THE INVERSE IN A
C  MR      IS THE MAXIMUM ROW DIMENSION OF A
C  MC      IS THE MAXIMUM COLUMN DIMENSION OF A
C  NR      IS THE NUMBER OF ROWS IN A
C  NC      IS THE NUMBER OF COLUMNS IN A
C  KZ      OUTPUT, KZ = 0 IMPLIES THAT A IS SINGULAR
C                  KZ = 1 IMPLIES THAT A IS NONSINGULAR
C  U       TEMPORARY WORKING STORAGE, DIMENSIONED MC BY MC
C  AFLAG   TEMPORARY WORKING STORAGE, DIMENSIONED MC
C  ATEMP   TEMPORARY WORKING STORAGE, DIMENSIONED MC
C
C-----------------------------------------------------------------------
      DIMENSION A(MR,MC),U(MC,MC),AFLAG(MC),ATEMP(MC)
      DOUBLE PRECISION A,U,AFLAG,ATEMP,FAC,DOT1,DOT2,TOL
     *,DOT,DSQRT
      KZ = 1
      TOL=.1D-21
      DO 10 I=1,NC
      DO 5 J=1,NC
    5 U(I,J)=0.0D0
   10 U(I,I)=1.0D0
      DO 12 L=1,NC
      FAC=DOT(MR,MC,NR,A,L,L)
      K=L
      IF(FAC.GT.TOL) GO TO 11
      DO 13 J=1,NR
   13 A(J,L)=0.0D0
      KZ=0
   12 AFLAG(L)=0.0D0
      GO TO 9997
   11 FAC = 1.0D0/DSQRT(FAC)
      L=K
      AFLAG(L)=1.0D0
      DO 15 I=1,NR
   15 A(I,L)=A(I,L)*FAC
      DO 20 I=1,NC
   20 U(I,L)=U(I,L)*FAC
C
      L1=L+1
      IF(L1.LE.NC) GO TO 21
      DO 22 I=1,NR
   22 A(I,L)=A(I,L)*FAC
      GO TO 9997
   21 DO 100 J=2,NC
      DOT1=DOT(MR,MC,NR,A,J,J)
      JM1=J-1
      DO 50 L=1,2
      DO 30 K=1,JM1
   30 ATEMP(K)=DOT(MR,MC,NR,A,J,K)
      DO 45 K=1,JM1
      DO 35 I=1,NR
   35 A(I,J)=A(I,J)-ATEMP(K)*A(I,K)*AFLAG(K)
      DO 40 I=1,NC
   40 U(I,J)=U(I,J)-ATEMP(K)*U(I,K)
   45 CONTINUE
   50 CONTINUE
      DOT2=DOT(MR,MC,NR,A,J,J)
      IF((DOT2/DOT1)-TOL) 55,55,70
   55 DO 60 I=1,JM1
      ATEMP(I)=0.00D0
      DO 60 K=1,I
   60 ATEMP(I)=ATEMP(I)+U(K,I)*U(K,J)
      DO 65 I=1,NR
      A(I,J)=0.0D0
      DO 65 K=1,JM1
   65 A(I,J)=A(I,J)-A(I,K)*ATEMP(K)*AFLAG(K)
      AFLAG(J)=0.0D0
      KZ = 0
C     FAC=DOT(MC,NC,U,J,J)
      FAC=DOT(MC,MC,NC,U,J,J)
      IF(FAC.GT.TOL) GO TO 66
      DO 67 I=1,NC
   67 U(I,J)=0.0D0
      GO TO 100
   66 FAC = 1.0D0/DSQRT(FAC)
      GO TO 75
   70 IF(DOT2.GT.TOL) GO TO 71
      KZ=0
      AFLAG(J)=0.0D0
      DO 72 I=1,NR
   72 A(I,J)=0.0D0
      GO TO 100
   71 AFLAG(J)=1.0D0
      FAC = 1.0D0/DSQRT(DOT2)
   75 DO 80 I=1,NR
   80 A(I,J)=A(I,J)*FAC
      DO 85 I=1,NC
   85 U(I,J)=U(I,J)*FAC
  100 CONTINUE
      DO 130 J=1,NC
      DO 130 I=1,NR
      FAC = 0.0D0
      DO 120 K=J,NC
  120 FAC=FAC+A(I,K)*U(J,K)
  130 A(I,J)=FAC
 9997 RETURN
      END
      REAL*8 FUNCTION DOT(MR,MC,NR,A,JC,KC)
      REAL*8 A(MR,MC)
C
C     COMPUTES THE INNER PRODUCT OF COLUMNS JC AND KC
C     OF MATRIX A.
      DOT = 0.0D0
      DO 5 I=1,NR
    5 DOT=DOT+A(I,JC)*A(I,KC)
      RETURN
      END

max_stars_repo_name: NASA-AMMOS/VICAR

max_issues_repo_name: NASA-AMMOS/VICAR

max_forks_repo_name: NASA-AMMOS/VICAR


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      SUBROUTINE DELMAT (N,NUC,NLC,M,A,B,IOPT)
!
!         THE SUBROUTINE DELMAT IS A DELFT-HYDRAULICS-LABORATORY
!         PRODUCT TO SOLVE SETS OF LINEAR EQUATIONS DESCRIBED BY
!         BAND MATRICES WITH LARGEST ELEMENTS ON THE DIAGONAL.
!         THE MAKER GIVES NO WARRANTY CONCERNING PROPER SOLUTIONS.
!         NOTHING OF THE CODE CAN BE REPRODUCED WITHOUT PERMISSION
!         OF THE MAKER
!
!         THE SUBROUTINE EXPECTS THE MATRIX "A" TO BE STORED IN A
!         ONE DIMENSIONAL EQUIVALENT OF THE TWO DIMENSIONAL MATRIX
!         REPRESENTATION. IN TWO DIMENSIONAL FORM THE STORAGE RULES
!         USED BY THE IMSL PACKAGE FOR BAND MATRICES MUST BE USED.
!         HOWEVER THE ROWS AND COLLUMNS MUST BE INTERCHANGED.
!         THE SAME HOLDS FOR THE FULL MATRIX OF KNOWN VECTORS "B".
!
!     SUBROUTINES CALLED  : SRSTOP, stops execution
!
!         NOTATION:                               CHANGED DURING
!                                                  CALCULATION
!
!      N    = ORDER OF MATRIX "A"                      NO
!      NUC  = NUMBER OF UPPER CODIAGONALS              NO
!      NLC  = NUMBER OF LOWER CODIAGONALS              NO
!      M    = NUMBER OF KNOWN VECTORS OF               NO
!             ORDER N IN "B"
!      A    = BANDMATRIX IN THE EQUIVALENT            YES
!             ONE DIMENSIONAL FORM
!      B    = MATRIX OF KNOWN VECTORS IN             YES/NO
!             ONE DIMENSIONAL FORM
!      IOPT = CALCULATION OPTION                       NO
!
!      OPTION 0 RETURNS THE LU-DECOMPOSITION IN MATRIX A, WHICH CAN
!                       BE USED AGAIN FOR NEW KNOWN VECTORS AND IT
!                       RETURNS THE SOLUTION FOR THE UNKNOWN VECTORS
!                       IN MATRIX B
!      OPTION 1 RETURNS ONLY THE DECOMPOSITION, B REMAINS UNCHANGED
!      OPTION 2 RETURNS THE SOLUTION FOR THE UNKNOWN VECTORS IN B
!                       IT NEEDS A PROPER DECOMPOSED AND STORED MATRIX
!                       A AS AN INPUT
!
      use timers

      DIMENSION A(*),B(*)
      integer(4) ithandl /0/
      if ( timon ) call timstrt ( "delmat", ithandl )
      NMUC = N - NUC
      NMLC = N - NLC
      NDM1 = NUC + NLC
      ND   = NDM1 + 1
      IF (IOPT.EQ.2) GOTO 1000
!
!           THE LU-DECOMPOSITION
!
!
!           K1 IS THE OUTER LOOP VARIABLE, COUNTING THE NUMBER OF
!              COLLUMNS WITH FULL LENGTH
!           L1 IS THE CORRESPONDING LIMIT VARIABLE
!           P  IS THE PIVOT ELEMENT. FOR THE ADVECTION DISPERSION
!              EQUATION THIS IS ALMOST ALWAYS THE LARGEST ELEMENT
!
      K1 = NLC + 1
      L1 = K1 + NMLC*ND
  100 P  = A(K1)
      IF ( ABS(P) .LT. 1.0E-35 ) THEN
          WRITE(6, '('' Matrix for DELMAT singular at element:'',I5)')
     *          K1/ND + 1
          CALL SRSTOP(1)
      ENDIF
!
!           K2 IS THE MIDDLE LOOP VARIABLE COUNTING THE NUMBER OF
!              ELEMENTS TO BE ELEMINATED
!           L2 IS THE CORRESPONDING LIMIT VARIABLE
!           F  IS THE MULTIPLICATION FACTOR FOR THE CORRECTION OF
!              THE ROW OF ELIMINATION. IT ENTERS A(K2) AS THE NEGATIVE
!              OF THE ELEMENT OF THE NEW LOWER TRIANGULAR MATRIX
!
      K2 = K1
      L2 = K1 + NLC*NDM1
  200 K2 = K2 + NDM1
      F  = A(K2)/P
      A(K2) = F
      K3 = K1
!
!           K4 IS THE INNER LOOP VARIABLE COUNTING THE ELEMENTS ON
!              THE ROW OF ELIMINATION
!           L4 IS THE CORRESPONDING LIMMIT VARIABLE
!           K3 IS THE CORRESPONDING COUNTER FOR THE ELEMENTS ON THE
!              ROW OF THE DIAGONAL ELEMENT
!
      K4 = K2
      L4 = K2 + NUC
  300 K3 = K3 + 1
      K4 = K4 + 1
      A(K4) = A(K4) - F * A(K3)
!
!           BOOKKEEPING OF THE LOOP VARIABLES
!
      IF (K4.LT.L4) GOTO 300
      IF (K2.LT.L2) GOTO 200
      K1 = K1 + ND
      IF (K1.LT.L1) GOTO 100
      IF (NLC.EQ.1) GOTO 700
!
!           THE COLLUMNS BECOME SHORTER, THE REST IS THE SAME
!
      N1 = NLC
      L1 = (N-1) * ND
  400 N1 = N1 - 1
      P  = A(K1)
      IF ( ABS(P) .LT. 1.0E-35 ) THEN
         WRITE(6,'('' Matrix for DELMAT singular at element:'',I5)')
     *         K1/ND + 1
         CALL SRSTOP(1)
      ENDIF
      K2 = K1
      L2 = K1 + N1*NDM1
  500 K2 = K2 + NDM1
      F  = A(K2)/P
      A(K2) = F
      K3 = K1
      K4 = K2
      L4 = K2 + NUC
  600 K3 = K3 + 1
      K4 = K4 + 1
      A(K4) = A(K4) - F * A(K3)
      IF (K4.LT.L4) GOTO 600
      IF (K2.LT.L2) GOTO 500
      K1 = K1 + ND
      IF (K1.LT.L1) GOTO 400
!
!           ENTRY FOR SUBSTITUTION OPTION
!
  700 IF(IOPT.EQ.1) goto 9999  !   RETURN
 1000 CONTINUE
!
!           THE FORWARD SUBSTITUTION HAS ESSENTIALLY THE SAME
!           STRUCTURE
!
!           K1 = OUTER LOOP COUNTER LOWER TRIANGULAR MATRIX
!           K2 = INNER LOOP COUNTER LOWER TRIANGULAR MATRIX
!           L1 AND L2 ARE THE LOOP LIMITS
!           K5 = OUTER LOOP COUNTER MATRIX OF KNOWN VECTORS
!           K4 = ROW COUNTER IN RELATION TO K5
!           K3 = INNER LOOP AND ROW COUNTER MATRIX OF KNOWN VECTORS
!           L3 = ROW LIMIT FOR ONE SUBSTITUTION ELEMENT "F"
!           F  = SUBSTITUTION ELEMENT OF LOWER TRIANGULAR MATRIX
!
      K1 = - NUC
      L1 = K1 + NMLC*ND
      K5 = - M
 1100 K1 = K1 + ND
      K5 = K5 + M
      K2 = K1
      L2 = K1 + NLC*NDM1
      K3 = K5 + M
      L3 = K3
 1200 K2 = K2 + NDM1
      F  = A(K2)
      L3 = L3 + M
      K4 = K5
 1300 K3 = K3 + 1
      K4 = K4 + 1
      B(K3) = B(K3) - F * B(K4)
      IF (K3.LT.L3) GOTO 1300
      IF (K2.LT.L2) GOTO 1200
      IF (K1.LT.L1) GOTO 1100
      IF (NLC.EQ.1) GOTO 2000
!
!           THE COLLUMNS BECOME SHORTER, THE REST IS THE SAME
!
      N1 = NLC
      L1 = (N-2) * ND
 1400 K1 = K1 + ND
      K5 = K5 + M
      K2 = K1
      N1 = N1 - 1
      L2 = K1 + N1*NDM1
      K3 = K5 + M
      L3 = K3
 1500 K2 = K2 + NDM1
      F  = A(K2)
      L3 = L3 + M
      K4 = K5
 1600 K3 = K3 + 1
      K4 = K4 + 1
      B(K3) = B(K3) - F * B(K4)
      IF (K3.LT.L3) GOTO 1600
      IF (K2.LT.L2) GOTO 1500
      IF (K1.LT.L1) GOTO 1400
!
!
!         BACKWARD SUBSTITUTION
!
 2000 K1 = N*ND + NLC + 1
      L1 = K1 - NMUC*ND
      K5 = N*M +1
 2100 K1 = K1 - ND
      L5 = K5 - M
      F  = A(K1)
 2200 K5 = K5 - 1
      B(K5) = B(K5)/F
      IF (K5.GT.L5) GOTO 2200
      K2 = K1
      L2 = K1 - NUC * NDM1
      K3 = K5
      L3 = K5
      K6 = K5 + M
 2300 K2 = K2 - NDM1
      F = A(K2)
      L3 = L3 - M
      K4 = K6
 2400 K3 = K3 - 1
      K4 = K4 - 1
      B(K3) = B(K3) - F * B(K4)
      IF (K3.GT.L3) GOTO 2400
      IF (K2.GT.L2) GOTO 2300
      IF (K1.GT.L1) GOTO 2100
      IF (NUC.EQ.1) GOTO 2850
      N1 = NUC
      L1 = 2*ND
 2500 K1 = K1 - ND
      L5 = K5 - M
      F  = A(K1)
      N1 = N1 - 1
 2600 K5 = K5 - 1
      B(K5) = B(K5)/F
      IF (K5.GT.L5) GOTO 2600
      K2 = K1
      L2 = K1 - N1 * NDM1
      K3 = K5
      L3 = K5
      K6 = K5 + M
 2700 K2 = K2 - NDM1
      F = A(K2)
      L3 = L3 - M
      K4 = K6
 2800 K3 = K3 - 1
      K4 = K4 - 1
      B(K3) = B(K3) - F * B(K4)
      IF (K3.GT.L3) GOTO 2800
      IF (K2.GT.L2) GOTO 2700
      IF (K1.GT.L1) GOTO 2500
 2850 K1 = K1 - ND
      L5 = K5 - M
      F  = A(K1)
 2900 K5 = K5 - 1
      B(K5) = B(K5)/F
      IF (K5.GT.L5) GOTO 2900

 9999 if ( timon ) call timstop ( ithandl )
      RETURN
      END

ext: f

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
module pffte
use types, only: dp
use ffte, only: dp_ffte, factor
use utils, only: stop_error
use mpi2, only: mpi_barrier
implicit none
private
public pfft3_init, pfft3, pifft3

interface
    subroutine PZFFT3DV(A,B,NX,NY,NZ,ICOMMY,ICOMMZ,NPUY,NPUZ,IOPT)
    import :: dp_ffte
    complex(dp_ffte), intent(inout) :: A(*), B(*)
    integer, intent(in) :: NX, NY, NZ, ICOMMY, ICOMMZ, NPUY, NPUZ, IOPT
    end subroutine
end interface

contains

subroutine pfft3_init(myid, comm_all, Ng, nsub)
integer, intent(in) :: myid, comm_all, Ng(:), nsub(:)
complex(dp) :: x(1), y(1)
integer :: commy, commz, ierr
integer :: Ng_local(3), F(3)
if (myid==0) then
    if (.not. all(nsub > 0)) then
        call stop_error("nsub must be positive")
    end if
    Ng_local = Ng / nsub
    if (nsub(1) /= 1) then
        call stop_error("nsub(1) must be equal to 1")
    end if
    if (.not. all(Ng_local * nsub == Ng)) then
        call stop_error("Ng must be equal to Ng_local * nsub")
    end if
    if (mod(Ng(1), nsub(2)) + mod(Ng(1), nsub(3)) + mod(Ng(2), nsub(3)) &
        + mod(Ng(3), nsub(2)) > 0) then
        call stop_error("Ng must be divisible by any permutation of nsub")
    end if
    call factor(Ng(1), F)
    if (Ng(1) /= 2**F(1)*3**F(2)*5**F(3)) then
        call stop_error("Ng(1) must be 2^a * 3^b * 5^c;  a,b,c=0,1,2,3,...")
    end if
    call factor(Ng(2), F)
    if (Ng(2) /= 2**F(1)*3**F(2)*5**F(3)) then
        call stop_error("Ng(2) must be 2^a * 3^b * 5^c;  a,b,c=0,1,2,3,...")
    end if
    call factor(Ng(3), F)
    if (Ng(3) /= 2**F(1)*3**F(2)*5**F(3)) then
        call stop_error("Ng(3) must be 2^a * 3^b * 5^c;  a,b,c=0,1,2,3,...")
    end if
end if
call mpi_barrier(comm_all, ierr)
! Only Ng is relevant if iopt == 0
! pzfft3dv divides by nsub inside, but does not store the result (the above
! checks make sure we do not divide by 0)
call pzfft3dv(x, y, Ng(1), Ng(2), Ng(3), commy, commz, nsub(2), nsub(3), 0)
end subroutine

subroutine pfft3(x, y, commy, commz, Ng, nsub)
complex(dp), intent(inout) :: x(:, :, :) ! input, will get destroyed
complex(dp), intent(out) :: y(:, :, :)   ! output
integer, intent(in) :: commy, commz ! communicators in y, z directions
integer, intent(in) :: Ng(:) ! Total (global) number of PW in each direction
integer, intent(in) :: nsub(:) ! Number of subdomains in each direction
call pzfft3dv(x, y, Ng(1), Ng(2), Ng(3), commy, commz, nsub(2), nsub(3), -1)
end subroutine

subroutine pifft3(x, y, commy, commz, Ng, nsub)
complex(dp), intent(inout) :: x(:, :, :) ! input, will get destroyed
complex(dp), intent(out) :: y(:, :, :)   ! output, will get divided by N
integer, intent(in) :: commy, commz ! communicators in y, z directions
integer, intent(in) :: Ng(:) ! Total (global) number of PW in each direction
integer, intent(in) :: nsub(:) ! Number of subdomains in each direction
call pzfft3dv(x, y, Ng(1), Ng(2), Ng(3), commy, commz, nsub(2), nsub(3), 1)
! pzfft3dv divides by N, so we multply by N here to cancel that
y = y * product(Ng)
end subroutine

end module

max_stars_repo_name: certik/hfsolver

max_issues_repo_name: certik/hfsolver

max_forks_repo_name: certik/hfsolver


! { dg-do compile }
! { dg-options "-fdump-tree-original" }
!
! PR fortran/41580
!
! Compile-time simplification of SAME_TYPE_AS and EXTENDS_TYPE_OF.

implicit none
type t1
  integer :: a
end type t1
type, extends(t1):: t11
  integer :: b
end type t11
type, extends(t11):: t111
  integer :: c
end type t111
type t2
  integer :: a
end type t2

type(t1) a1
type(t11) a11
type(t2) a2
class(t1), allocatable :: b1
class(t11), allocatable :: b11
class(t2), allocatable :: b2

logical, parameter :: p1 = same_type_as(a1,a2)  ! F
logical, parameter :: p2 = same_type_as(a2,a1)  ! F
logical, parameter :: p3 = same_type_as(a1,a11) ! F
logical, parameter :: p4 = same_type_as(a11,a1) ! F
logical, parameter :: p5 = same_type_as(a11,a11)! T
logical, parameter :: p6 = same_type_as(a1,a1)  ! T

if (p1 .or. p2 .or. p3 .or. p4 .or. .not. p5 .or. .not. p6) call should_not_exist()

if (same_type_as(b1,b1)   .neqv. .true.) call should_not_exist()

! Not (trivially) compile-time simplifiable:
if (same_type_as(b1,a1)  .neqv. .true.) call abort()
if (same_type_as(b1,a11) .neqv. .false.) call abort()
allocate(t1 :: b1)
if (same_type_as(b1,a1)  .neqv. .true.) call abort()
if (same_type_as(b1,a11) .neqv. .false.) call abort()
deallocate(b1)
allocate(t11 :: b1)
if (same_type_as(b1,a1)  .neqv. .false.) call abort()
if (same_type_as(b1,a11) .neqv. .true.) call abort()
deallocate(b1)

! .true. -> same type
if (extends_type_of(a1,a1)   .neqv. .true.) call should_not_exist()
if (extends_type_of(a11,a11) .neqv. .true.) call should_not_exist()
if (extends_type_of(a2,a2)   .neqv. .true.) call should_not_exist()

! .false. -> type compatibility possible
if (extends_type_of(a1,a2)  .neqv. .false.) call should_not_exist()
if (extends_type_of(a2,a1)  .neqv. .false.) call should_not_exist()
if (extends_type_of(a11,a2) .neqv. .false.) call should_not_exist()
if (extends_type_of(a2,a11) .neqv. .false.) call should_not_exist()

if (extends_type_of(b1,b2)  .neqv. .false.) call should_not_exist()
if (extends_type_of(b2,b1)  .neqv. .false.) call should_not_exist()
if (extends_type_of(b11,b2) .neqv. .false.) call should_not_exist()
if (extends_type_of(b2,b11) .neqv. .false.) call should_not_exist()

if (extends_type_of(b1,a2)  .neqv. .false.) call should_not_exist()
if (extends_type_of(b2,a1)  .neqv. .false.) call should_not_exist()
if (extends_type_of(b11,a2) .neqv. .false.) call should_not_exist()
if (extends_type_of(b2,a11) .neqv. .false.) call should_not_exist()

if (extends_type_of(a1,b2)  .neqv. .false.) call should_not_exist()
if (extends_type_of(a2,b1)  .neqv. .false.) call should_not_exist()
if (extends_type_of(a11,b2) .neqv. .false.) call should_not_exist()
if (extends_type_of(a2,b11) .neqv. .false.) call should_not_exist()

! type extension possible, compile-time checkable
if (extends_type_of(a1,a11) .neqv. .false.) call should_not_exist()
if (extends_type_of(a11,a1) .neqv. .true.) call should_not_exist()

if (extends_type_of(b1,a1)   .neqv. .true.) call should_not_exist()
if (extends_type_of(b11,a1)  .neqv. .true.) call should_not_exist()
if (extends_type_of(b11,a11) .neqv. .true.) call should_not_exist()

if (extends_type_of(a1,b11)  .neqv. .false.) call should_not_exist()

! Special case, simplified at tree folding:
if (extends_type_of(b1,b1)   .neqv. .true.) call abort()

! All other possibilities are not compile-time checkable
if (extends_type_of(b11,b1)  .neqv. .true.) call abort()
if (extends_type_of(b1,b11)  .neqv. .false.) call abort()
if (extends_type_of(a11,b11) .neqv. .true.) call abort()

allocate(t11 :: b11)
if (extends_type_of(a11,b11) .neqv. .true.) call abort()
deallocate(b11)

allocate(t111 :: b11)
if (extends_type_of(a11,b11) .neqv. .false.) call abort()
deallocate(b11)

allocate(t11 :: b1)
if (extends_type_of(a11,b1) .neqv. .true.) call abort()
deallocate(b1)

allocate(t11::b1)
if (extends_type_of(b1,a11) .neqv. .true.) call abort()
deallocate(b1)

allocate(b1,source=a11)
if (extends_type_of(b1,a11) .neqv. .true.) call abort()
deallocate(b1)

allocate( b1,source=a1)
if (extends_type_of(b1,a11) .neqv. .false.) call abort()
deallocate(b1)

end

! { dg-final { scan-tree-dump-times "abort" 16 "original" } }
! { dg-final { scan-tree-dump-times "should_not_exist" 0 "original" } }


text:
module model_parser_module

  ! read in an initial model and return arrays with the model data.
  ! take care to match the species available in the model file to
  ! those defined by the network

  ! the model file is assumed to be of the follow form:
  ! # npts = 896
  ! # num of variables = 6
  ! # density
  ! # temperature
  ! # pressure
  ! # carbon-12
  ! # oxygen-16
  ! # magnesium-24
  ! 195312.5000  5437711139.  8805500.952   .4695704813E+28  0.3  0.7  0
  ! 585937.5000  5410152416.  8816689.836  0.4663923963E+28  0.3  0.7  0
  !
  ! we read in the number of variables and their order and use this to
  ! map them into the model_state array.  We ignore anything other than
  ! density, temperature, pressure and composition.
  !
  ! composition is assumed to be in terms of mass fractions

  use parallel, only: parallel_IOProcessor
  use network
  use bl_types

  implicit none

  ! integer keys for indexing the model_state array
  integer, parameter :: nvars_model = 3 + nspec
  integer, parameter :: idens_model = 1
  integer, parameter :: itemp_model = 2
  integer, parameter :: ipres_model = 3
  integer, parameter :: ispec_model = 4

  ! number of points in the model file
  integer, save :: npts_model

  ! arrays for storing the model data
  real (kind=dp_t), allocatable, save :: model_state(:,:)
  real (kind=dp_t), allocatable, save :: model_r(:)

  ! model_initialized will be .true. once the model is read in and the
  ! model data arrays are initialized and filled
  logical, save :: model_initialized = .false.

  integer, parameter :: MAX_VARNAME_LENGTH=80

  public :: read_model_file, close_model_file

contains

  subroutine read_model_file(model_file)

    use bl_constants_module
    use bl_error_module

    character(len=*), intent(in   ) :: model_file

    ! local variables
    integer :: nvars_model_file
    integer :: ierr

    integer :: i, j, comp

    real(kind=dp_t), allocatable :: vars_stored(:)
    character(len=MAX_VARNAME_LENGTH), allocatable :: varnames_stored(:)
    logical :: found_model, found_dens, found_temp, found_pres
    logical :: found_spec(nspec)
    integer :: ipos
    character (len=256) :: header_line

    ! open the model file
    open(99,file=trim(model_file),status='old',iostat=ierr)

    if (ierr .ne. 0) then
       print *,'Couldnt open model_file: ',model_file
       call bl_error('Aborting now -- please supply model_file')
    end if

    ! the first line has the number of points in the model
    read (99, '(a256)') header_line
    ipos = index(header_line, '=') + 1
    read (header_line(ipos:),*) npts_model

    ! now read in the number of variables
    read (99, '(a256)') header_line
    ipos = index(header_line, '=') + 1
    read (header_line(ipos:),*) nvars_model_file

    allocate (vars_stored(nvars_model_file))
    allocate (varnames_stored(nvars_model_file))

    ! now read in the names of the variables
    do i = 1, nvars_model_file
       read (99, '(a256)') header_line
       ipos = index(header_line, '#') + 1
       varnames_stored(i) = trim(adjustl(header_line(ipos:)))
    enddo

    ! allocate storage for the model data
    allocate (model_state(npts_model, nvars_model))
    allocate (model_r(npts_model))

887 format(78('-'))
889 format(a60)

    if ( parallel_IOProcessor()) then
       write (*,889) ' '
       write (*,887)
       write (*,*)   'reading initial model'
       write (*,*)   npts_model, 'points found in the initial model file'
       write (*,*)   nvars_model_file, ' variables found in the initial model file'
    endif

    ! start reading in the data
    do i = 1, npts_model
       read(99,*) model_r(i), (vars_stored(j), j = 1, nvars_model_file)

       model_state(i,:) = ZERO

       ! make sure that each of the variables that MAESTRO cares about
       ! are found
       found_dens = .false.
       found_temp = .false.
       found_pres = .false.
       found_spec(:) = .false.

       do j = 1,nvars_model_file

          ! keep track of whether the current variable from the model
          ! file is one that MAESTRO cares about
          found_model = .false.

          if (varnames_stored(j) == "density") then
             model_state(i,idens_model) = vars_stored(j)
             found_model = .true.
             found_dens  = .true.

          else if (varnames_stored(j) == "temperature") then
             model_state(i,itemp_model) = vars_stored(j)
             found_model = .true.
             found_temp  = .true.

          else if (varnames_stored(j) == "pressure") then
             model_state(i,ipres_model) = vars_stored(j)
             found_model = .true.
             found_pres  = .true.

          else
             do comp = 1, nspec
                if (varnames_stored(j) == spec_names(comp)) then
                   model_state(i,ispec_model-1+comp) = vars_stored(j)
                   found_model = .true.
                   found_spec(comp) = .true.
                   exit
                endif
             enddo
          endif

          ! is the current variable from the model file one that we
          ! care about?
          if (.NOT. found_model .and. i == 1) then
             if ( parallel_IOProcessor() ) then
                print *, 'WARNING: variable not found: ', &
                     trim(varnames_stored(j))
             end if
          endif

       enddo   ! end loop over nvars_model_file

       ! were all the variables we care about provided?
       if (i == 1) then
          if (.not. found_dens) then
             if ( parallel_IOProcessor() ) then
                print *, 'WARNING: density not provided in inputs file'
             end if
          endif

          if (.not. found_temp) then
             if ( parallel_IOProcessor() ) then
                print *, 'WARNING: temperature not provided in inputs file'
             end if
          endif

          if (.not. found_pres) then
             if ( parallel_IOProcessor() ) then
                print *, 'WARNING: pressure not provided in inputs file'
             end if
          endif

          do comp = 1, nspec
             if (.not. found_spec(comp)) then
                if ( parallel_IOProcessor() ) then
                   print *, 'WARNING: ', trim(spec_names(comp)), &
                        ' not provided in inputs file'
                end if
             endif
          enddo
       endif

    end do   ! end loop over npts_model

    close(99)

    model_initialized = .true.

    deallocate(vars_stored,varnames_stored)

  end subroutine read_model_file

  function get_model_npts(model_file)

    integer :: get_model_npts

    ! look in the model file and return the number of points
    character(len=256), intent(in   ) :: model_file

    character (len=256) :: header_line
    integer :: ipos

    open(99,file=model_file)

    ! the first line has the number of points in the model
    read (99, '(a256)') header_line
    ipos = index(header_line, '=') + 1
    read (header_line(ipos:),*) get_model_npts

    close(99)

    return

  end function get_model_npts

  subroutine close_model_file

    if (model_initialized) then
       deallocate(model_r)
       deallocate(model_state)
       npts_model = -1
       model_initialized = .false.
    endif
  end subroutine close_model_file

end module model_parser_module

max_stars_repo_name: yingtchen/Castro

max_issues_repo_name: yingtchen/Castro

max_forks_repo_name: yingtchen/Castro


text:
!
! This file is released under terms of BSD license
! See LICENSE file for more information
!
program openacc_continuation
  implicit none

  integer :: i
  real, dimension(100)  :: a
  real, dimension(100)  :: bvariable, cvariable, dvariable, evariable, &
                           fvariable, gvariable, hvariable, ivariable, &
                           jvariable, kvariable, lvariable

  !$claw acc data copyin(a) &
  !$claw copyin(bvariable, cvariable, dvariable, evariable, fvariable) &
  !$claw copyin(gvariable, hvariable, ivariable, jvariable) &
  !$claw copyin(kvariable, lvariable)
  !$claw acc parallel
  !$claw acc loop gang vector
  do i=1, 100
    a(i) = bvariable(i)*2.0
  end do
  !$claw acc end parallel
  !$claw acc end data
end program openacc_continuation

max_stars_repo_name: FrostyMike/claw-compiler

max_issues_repo_name: FrostyMike/claw-compiler

max_forks_repo_name: FrostyMike/claw-compiler


text:
Describe Users/KatherineLucas here.

20110809 17:20:15 nbsp Hi Katherine, and Welcome to the Wiki.  Please take a minute to read over Welcome to the Wiki/Business Owner this pageit highlights some important points about the wiki.  First and foremost, as a nonprofit, its prohibited from the IRS by being used for promotion of forprofit businesses.  Your edits to Hyatt Place UC Davis read just like a flyer or website promoting the hotelwhat we call glossy on the Wiki.  That sort of thing is specifically prohibited, so I undid or altered a number of the changes you made to be less promotional.

Thanks for editing, and feel free to ask anywhere if you have any questions. Users/TomGarberson

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search