$NOFLOATCALLS
      SUBROUTINE STBCVN
      IMPLICIT REAL*8(A-H,O-Z)
C----------------shared variables------------
      REAL*8 B,TBCVNV,FVAL,TBFV,DET,FBASE,FIT
      INTEGER NELTS,NOBS,MAXELT,MAXOBS,MOBSYM
      LOGICAL LNFLAG
      PARAMETER(MAXELT=15,MAXOBS=120,MOBSYM=MAXOBS*(MAXOBS+1)/2)
      COMMON /WORK/ DET, B(MAXELT,MAXOBS),TBCVNV(MOBSYM),
     1       FBASE,FVAL(MAXOBS),TBFV(MAXOBS),NELTS,NOBS,FIT,LNFLAG
C---------------------------------------------
C     This routine takes B as given, computes TBCVNV from it.
C     TBCVNV is inverse of est.'d covariance matrix
C     across data points.
C---------------------------------------------
      REAL*8 DIST, X, DOT
      LOGICAL NOTOK,ANSYES
      INTEGER NZEROS, II, I, J, LUN
      DATA PRECIS /1.D-15/
      NZEROS=0
	II=0
	DO 10 I=1,NOBS
	    DO 11 J=1,I-1
		II=II+1
11	   	TBCVNV(II)=DIST(B(1,J),B(1,I))
		II=II+1
	    TBCVNV(II)=1.0
10	CONTINUE
      CALL CHOLES(TBCVNV,NOBS,NOTOK)
	IF(NOTOK) THEN
	    WRITE(*,*) 'Accuracy problems in setup. Likelihood dubious.'
	ENDIF
      DET=0.0
      II=0
      DO 15 I=1,NOBS
  		II=II+I
  		IF(TBCVNV(II).GE.PRECIS) THEN
            DET=DET+LOG(TBCVNV(II))		  
         ELSE
            NZEROS = NZEROS+1
         ENDIF
15    CONTINUE
C----------------------------------------------------
C     Note that DET is log det of inverse of scaled VCV matrix
C-----------------------------------------------------
      DET=2.*DET
      CALL TRIINV(TBCVNV,NOBS)
      CALL TRISQ(TBCVNV,NOBS)
      CALL SYMVEC(TBCVNV,FVAL,TBFV,NOBS)
C      -------------------------------------------------
C      FIT is -2 times log of Lhood concentrated w.r.t. a scale
C      parameter on the covariance matrix.
C      --------------------------------------------------
      X=DOT(TBFV,FVAL,NOBS)
      IF (X .GT.1.D-80) THEN
         FIT=-(NOBS-NZEROS)*LOG(DOT(TBFV,FVAL,NOBS)/(NOBS-NZEROS))-DET
      ELSE
         FIT = 9.9999999999D99
      ENDIF
	WRITE(*,*) 'Write out projection matrix?'
	IF(ANSYES()) THEN
	    WRITE(*,*) 'To screen?'
	    IF(ANSYES()) THEN
		LUN=0
	    ELSE
		CALL DEFOUT(6)
		LUN=6
	    ENDIF
	    CALL SYMWRI(TBCVNV,LUN,NOBS)
	    IF(LUN.EQ.6) CLOSE(6)
	ENDIF
      RETURN
      END
  
      FUNCTION SMTH(A)
      IMPLICIT REAL*8(A-H,O-Z)
C----------------shared variables------------
      REAL*8 B,TBCVNV,FVAL,TBFV, DET, FBASE, FIT
      INTEGER NELTS,NOBS,MAXELT,MAXOBS,MOBSYM
      LOGICAL LNFLAG
      PARAMETER(MAXELT=15,MAXOBS=120,MOBSYM=MAXOBS*(MAXOBS+1)/2)
      COMMON /WORK/ DET, B(MAXELT,MAXOBS),TBCVNV(MOBSYM),
     1       FBASE,FVAL(MAXOBS),TBFV(MAXOBS),NELTS,NOBS,FIT,LNFLAG
C---------------------------------------------
      REAL*8 WORVEC(MAXOBS),CORVEC(MAXOBS),A(*), DIST, SMTH, DOT
      INTEGER I
      DO 10 I=1,NOBS
10       CORVEC(I)= DIST(A,B(1,I))
      CALL SYMVEC(TBCVNV,CORVEC,WORVEC,NOBS)
      SMTH=DOT(WORVEC,FVAL,NOBS)+FBASE
      RETURN
      END

      SUBROUTINE PTNVUP(K,B,DET,N,OK)
      IMPLICIT REAL*8(A-H,O-Z)
C------------------------------------------------
C    |K   0|                                       |A  B| -1
C    |     |  = Inverse of lower triang. sqrt. of  |    |  .  On entry, K
C    |L'  M|                                       |B' C|
C
C    is the transpose of the inverse of the lower triangular square root of
C    A. )  The B argument is B in the formula,
C    augmented by C at the end.  On return K is an order N+1 matrix.
C    This is set up to use the unpivoted version of TRIVEC and VECTRI.
C-------------------------------------------------
      PARAMETER(PRECIS=1.D-30)
      INTEGER BLOC
      LOGICAL OK
      REAL*8 K(*),B(*),DET, C, DOT
      BLOC=N*(N+1)/2+1
      C=B(N+1)
      CALL TRIVEC(K,B,K(BLOC),N)
      C = C-DOT(K(BLOC),K(BLOC),N)
      IF (C.GT.PRECIS) THEN
         K(BLOC+N)=1./SQRT(C)
         CALL VECTRI(K(BLOC),K,B,N)
         CALL MATMOV(B,K(BLOC),N)
         CALL SCALAR(K(BLOC),-K(BLOC+N),N)
         DET=DET*K(BLOC+N)
         OK=.TRUE.
      ELSE
         OK=.FALSE.
         IF (K(BLOC+N).LT.-PRECIS) THEN
            WRITE(*,'('' NEGATIVE DEFINITE IN PTNVUP'')')
         ELSE
             WRITE(*,'('' PRECISION PROBLEM IN PTNVUP'')')
         ENDIF
      ENDIF
      RETURN
      END


