C @(#)prd.for	19.1 (ES0-DMD) 02/25/03 14:29:02
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
CC
CC  *RANGS* : SORTS A VECTOR OF LENGTH 'NCAS'.
CC
      SUBROUTINE RANGS(AM,LAME)

C       INCLUDE 'MID_REL_INCL:implicit.inc'
C      DIMENSION JLV(30),JRV(30),AM(LAME)

      INTEGER J,JNC,JNDL,JR,JSS,JTWE,LAME
      REAL AMM,XX

      INTEGER JLV(30),JRV(30)
      REAL AM(LAME)

      JSS=1
      JLV(1)=1
      JRV(1)=LAME
  10  JNDL=JLV(JSS)
      JR=JRV(JSS)
      JSS=JSS-1
  20  JNC=JNDL
      J=JR
      JTWE=(JNDL+JR)/2
      XX=AM(JTWE)
  30  IF(AM(JNC).GE.XX)GOTO 40
      JNC=JNC+1
      GOTO 30
  40  IF(XX.GE.AM(J))GOTO 50
      J=J-1
      GOTO 40
  50  IF(JNC.GT.J)GOTO 60
      AMM=AM(JNC)
      AM(JNC)=AM(J)
      AM(J)=AMM
      JNC=JNC+1
      J=J-1
  60  IF (JNC.LE.J) GOTO 30
      IF ((J-JNDL).LT.(JR-JNC)) GOTO 80
      IF (JNDL.GE.J) GOTO 70
      JSS=JSS+1
      JLV(JSS)=JNDL
      JRV(JSS)=J
  70  JNDL=JNC
      GOTO 100
  80  IF (JNC.GE.JR) GOTO 90
      JSS=JSS+1
      JLV(JSS)=JNC
      JRV(JSS)=JR
  90  JR=J
  100 IF (JNDL.LT.JR) GOTO 20
      IF (JSS.NE.0) GOTO 10
      RETURN
      END
CC
CC  *RDUAL* : CALCULATES THE RESIDUALS OF ALL CASES
CC
       SUBROUTINE RDUAL(AA,JKEUS,JDA,NCAS,NVAR,JCST,JPRT,NVAD,LUB,PREC,
     1 JREG,X,Y,RESDU,WEIGHTS,XMED,XMAD,NZWE,AVW,JPLT,AW,NMVAL,
     1 NDXX,NDXY,JHEAD,LAB)

C       INCLUDE 'MID_REL_INCL:implicit.inc'

C       DIMENSION AA(JDA),X(NDXX,NDXY),Y(NDXY),RESDU(NDXY),WEIGHTS(NDXY)
C       DIMENSION XMED(NDXX),XMAD(NDXX),AW(NDXY),NMVAL(NDXY)
       INTEGER J,JCST,JDA,JKEUS,JNC,JPL2,JPLT,JPRT,JREG
       INTEGER NCAS,NDXX,NDXY,NVAD,NVAR,NZWE,LUB
       REAL AL,AVW,BBB,PREC,YJ
 
       REAL AA(JDA),X(NDXX,NDXY),Y(NDXY),RESDU(NDXY),WEIGHTS(NDXY)
       REAL XMED(NDXX),XMAD(NDXX),AW(NDXY)
       INTEGER NMVAL(NDXY)


       CHARACTER*10 LAB(NDXX)
       CHARACTER*60 JHEAD
       JPL2=1
       AL=NCAS
       IF(JPRT.EQ.0) GOTO 10
       IF (NVAR.EQ.1.AND.JCST.EQ.1) THEN
          IF (JKEUS.NE.2) WRITE(LUB,8000) LAB(NVAD)
          IF (JKEUS.EQ.2) WRITE(LUB,8010) LAB(NVAD)
                                    ELSE
          IF (JKEUS.NE.2) WRITE(LUB,8020) LAB(NVAD),LAB(NVAD)
          IF (JKEUS.EQ.2) WRITE(LUB,8030) LAB(NVAD),LAB(NVAD)
       ENDIF
 10    DO 20 JNC=1,NCAS
       IF (.NOT.(NVAR.EQ.1.AND.JCST.EQ.1)) GOTO 40
       BBB=Y(JNC)-AA(1)
       GOTO 60
 40    AW(JNC)=0.0
       DO 30 J=1,NVAR
 30    AW(JNC)=AW(JNC)+AA(J)*(X(J,JNC)*XMAD(J)+XMED(J))
       YJ=Y(JNC)*XMAD(NVAD)+XMED(NVAD)
       BBB=YJ-AW(JNC)
 60    IF (AA(NVAD).GT.PREC) GOTO 70
       IF (ABS(BBB).LT.PREC) BBB=0.0
       RESDU(JNC)=BBB
       GOTO 50
 70    RESDU(JNC)=BBB/AA(NVAD)
 50    IF (JPRT.EQ.0) GOTO 20
       IF (NVAR.EQ.1.AND.JCST.EQ.1) GOTO 80
       IF (AA(NVAD).GT.PREC.AND.JKEUS.NE.2)
     1 WRITE(LUB,8040) YJ,AW(JNC),BBB,NMVAL(JNC),RESDU(JNC)
       IF (AA(NVAD).LE.PREC)
     1 WRITE(LUB,8050) YJ,AW(JNC),BBB,NMVAL(JNC)
       IF (AA(NVAD).GT.PREC.AND.JKEUS.EQ.2)
     1 WRITE(LUB,8040) YJ,AW(JNC),BBB,NMVAL(JNC),RESDU(JNC),
     1 WEIGHTS(JNC)
       GOTO 20
 80    IF (JKEUS.NE.2) WRITE(LUB,8060) Y(JNC),BBB,JNC,RESDU(JNC)
       IF (JKEUS.EQ.2) WRITE(LUB,8060) Y(JNC),BBB,JNC,RESDU(JNC),
     1 WEIGHTS(JNC)
 20    CONTINUE
       IF(JPLT.LT.1) GOTO 90
       IF (AA(NVAD).LE.PREC) JPL2=0
       CALL GRAF(AW,RESDU,NCAS,0,JPL2,JPLT,NMVAL,NDXY,LUB,JREG,
     1 JHEAD,NDXX,NVAD,LAB)
       IF (JPLT.LT.3) GOTO 90
       JPLT=2
       CALL GRAF(AW,RESDU,NCAS,0,JPL2,JPLT,NMVAL,NDXY,LUB,JREG,
     1 JHEAD,NDXX,NVAD,LAB)
       JPLT=3
 90    CONTINUE
 8000  FORMAT(//11X,'OBSERVED',12X,'RESIDUAL',
     1 '   NO','  RES/SC'/9X,A10/)
 8010  FORMAT(//11X,'OBSERVED',12X,'RESIDUAL',
     1 '   NO','  RES/SC',' WEIGHT'/9X,A10/)
 8020  FORMAT(//11X,'OBSERVED',11X,'ESTIMATED',12X,'RESIDUAL',
     1 3X,'NO','  RES/SC'/9X,A10,10X,A10/)
 8030  FORMAT(//11X,'OBSERVED',11X,'ESTIMATED',12X,'RESIDUAL',
     1 3X,'NO','  RES/SC',' WEIGHT'/9X,A10,10X,A10/)
 8040  FORMAT(1X,F18.5,2(1X,F19.5),I5,1X,F7.2,F6.1)
 8050  FORMAT(1X,F18.5,2(1X,F19.5),I5,'   *****')
 8060  FORMAT(1X,F18.5,1X,F19.5,I5,1X,F7.2,F6.1)
       RETURN
       END
CC
CC  *SHHLF* : CALCULATES THE LMS IN THE ONE-DIMENSIONAL CASE.
CC
       SUBROUTINE SHHLF(W,NCAS,JQU,SLUTN,BSTD,PREC)

C       INCLUDE 'MID_REL_INCL:implicit.inc'
C       DIMENSION W(NCAS)
       INTEGER JNC,JPK,JQU,JTEL,K,LNK,NCAS
       REAL ABSW,BSTD,PREC,SLUTN,TEL,VERSC

       REAL W(NCAS)

       K=JQU-1
       VERSC=W(K+1)-W(1)
       SLUTN=(W(1)+W(K+1))/2.0
       BSTD=((W(K+1)-W(1))/2.0)
       IF (NCAS.EQ.2) RETURN
       JTEL=1
       LNK=NCAS-K
       DO 10 JNC=2,LNK
       JPK=JNC+K
       ABSW=W(JPK)-W(JNC)
       IF ( ABS(ABSW-VERSC).GT.PREC) GOTO 20
       JTEL=JTEL+1
       SLUTN=SLUTN+(W(JNC)+W(JPK))/2.0
       GOTO 10
 20    IF (ABSW.LT.VERSC) THEN
       VERSC=ABSW
       JTEL=1
       BSTD=((W(JPK)-W(JNC))/2.0)
       SLUTN=(W(JNC)+W(JPK))/2.0
       ENDIF
 10    CONTINUE
       TEL=JTEL
       SLUTN=SLUTN/TEL
       RETURN
       END
CC
CC  *TRC* : RETRANSFORMATION OF THE VARIANCE-COVARIANCE MATRIX
CC
       SUBROUTINE TRC(H,DA,NMXV,NDXX,NDXY,NVAR,JCST,NVSB,NVAD,
     1 XMED,XMAD)

C       INCLUDE 'MID_REL_INCL:implicit.inc'
C       DIMENSION H(NMXV,NDXX),DA(NMXV)
C       DOUBLE PRECISION H,XMP2,HNN
C       DIMENSION XMED(NDXX),XMAD(NDXX)
       INTEGER  J,JCST,JU,K,K2,NDXX,NDXY,NMXV,NVAD,NVAR,NVSB

       REAL DA(NMXV)
       DOUBLE PRECISION H(NMXV,NDXX),XMP2,HNN
       REAL XMED(NDXX),XMAD(NDXX)

       XMP2=DBLE(XMAD(NVAD))*DBLE(XMAD(NVAD))
       IF (JCST.EQ.0) THEN
          DO 10 J=1,NVAR
          DO 20 K=1,J
          H(J,K)=H(J,K)*(XMP2/(DBLE(XMAD(J))*DBLE(XMAD(K))))
 20          CONTINUE
          DA(J)=DSQRT(H(J,J))
 10          CONTINUE
       ELSE
          DO 25 J=1,NVAR
          H(J,NVAD)=H(J,J)
 25          CONTINUE
          DO 30 J=1,NVAR
          DO 40 K=1,J
          H(J,K)=H(J,K)*XMP2/(DBLE(XMAD(J))*DBLE(XMAD(K)))
 40          CONTINUE
          DA(J)=DSQRT(H(J,J))
 30          CONTINUE
          DO 50 K=1,NVSB
          H(NVAR,K)=H(K,NVAR)*XMP2/DBLE(XMAD(K))
          DO 60 K2=1,NVAR
          IF (K.EQ.K2) THEN
            H(NVAR,K)=H(NVAR,K)-DBLE(XMED(K))*XMP2/
     1                      (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K2,NVAD)
            GOTO 60
          ENDIF
          IF (K.LT.K2) THEN
               H(NVAR,K)=H(NVAR,K)-(DBLE(XMED(K2))*XMP2)/
     1                      (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K,K2)
                        ELSE
               H(NVAR,K)=H(NVAR,K)-DBLE(XMED(K2))*XMP2/
     1                      (DBLE(XMAD(K2))*DBLE(XMAD(K)))*H(K2,K)
          ENDIF
 60          CONTINUE
 50          CONTINUE
          H(NVAR,NVAR)=H(NVAR,NVAD)*XMP2
          DO 70 K=1,NVAR
          H(NVAR,NVAR)=H(NVAR,NVAR)+
     1            (DBLE(XMED(K))*DBLE(XMED(K)))*XMP2/
     1                         (DBLE(XMAD(K))*DBLE(XMAD(K)))*H(K,NVAD)
 70          CONTINUE
          DO 80 K=1,NVAR
          IF (K.NE.NVAR) THEN
             H(NVAR,NVAR)=H(NVAR,NVAR)-2.0D0*XMP2*DBLE(XMED(K))/
     1                                    (DBLE(XMAD(K)))*H(K,NVAR)
                         ELSE
             H(NVAR,NVAR)=H(NVAR,NVAR)-2.0D0*XMP2*DBLE(XMED(K))/
     1                                    (DBLE(XMAD(K)))*H(NVAR,NVAD)
          ENDIF
 80          CONTINUE
          DO 90 J=1,NVSB
          JU=J+1
          DO 90 K=JU,NVAR
          HNN=2.0D0*DBLE(XMED(J))*DBLE(XMED(K))*XMP2
          H(NVAR,NVAR)=H(NVAR,NVAR)+HNN/
     1                              (DBLE(XMAD(J))*DBLE(XMAD(K)))*H(J,K)
 90          CONTINUE
          DA(NVAR)=DSQRT(H(NVAR,NVAR))
       ENDIF
       RETURN
       END
CC
CC  *SCHCV* : PRINTS THE VARIANCE-COVARIANCE MATRIX OF THE
CC            REGRESSION COEFFICIENTS.
CC
       SUBROUTINE SCHCV(NVAR,NDXX,NMXV,H,LUB)

C       INCLUDE 'MID_REL_INCL:implicit.inc'
       INTEGER J,K,NDXX,NMXV,NVAR,LUB

       DOUBLE PRECISION H(NMXV,NDXX)
       WRITE(LUB,9000)
       DO 10 J=1,NVAR
       WRITE(LUB,9010) (H(J,K),K=1,J)
 10    CONTINUE
 9000  FORMAT(//' VARIANCE - COVARIANCE MATRIX = '/)
 9010  FORMAT(5(5X,D10.4))
       RETURN
       END
CC
CC  *RTRAN* : RETRANSFORMATION OF THE REGRESSION COEFFICIENTS.
CC
       SUBROUTINE RTRAN(NVAR,JCST,NVSB,NVAD,NDXX,XMED,XMAD,AA,JAL,FCKW)

C       INCLUDE 'MID_REL_INCL:implicit.inc'
C       DIMENSION AA(JAL),XMED(NDXX),XMAD(NDXX)

       INTEGER J,JAL,JCST,NDXX,NVAD,NVAR,NVSB
       REAL  FCKW

       REAL AA(JAL),XMED(NDXX),XMAD(NDXX)

       IF (NVAR.LE.1) THEN
         AA(1)=AA(1)*XMAD(NVAD)/XMAD(1)
         GOTO 30
       ENDIF
       DO 10 J=1,NVSB
 10    AA(J)=AA(J)*XMAD(NVAD)/XMAD(J)
       IF (JCST.EQ.0) THEN
         AA(NVAR)=AA(NVAR)*XMAD(NVAD)/XMAD(NVAR)
                      ELSE
         AA(NVAR)=AA(NVAR)*XMAD(NVAD)
         DO 20 J=1,NVSB
 20         AA(NVAR)=AA(NVAR)-AA(J)*XMED(J)
         AA(NVAR)=AA(NVAR)+XMED(NVAD)
       ENDIF
 30    FCKW=FCKW*(XMAD(NVAD)*XMAD(NVAD))
       RETURN
       END
CC
CC  *RANPN* : RANDOM NUMBER GENERATOR .
CC
       SUBROUTINE RANPN(NCAS,NVAR,JNDVC,NDXX,JEY,JLP,MAXRP)
C       INCLUDE 'MID_REL_INCL:implicit.inc'
C       DIMENSION JNDVC(NDXX)
       INTEGER J,JEY,JEYQT,JLP,JNC,JNCM,JRAN,MAXRP,NCAS,NDXX,NVAR
       REAL AL,RNULE,RY

       INTEGER JNDVC(NDXX)

CC     WE PROGRAMMED THIS RANDOM NUMBER GENERATOR OURSELVES BECAUSE
CC     WE WANTED IT TO BE MACHINE INDEPENDENT. IT SHOULD RUN ON MOST
CC     COMPUTERS BECAUSE THE LARGEST INTEGER USED IS LESS THAN 2**30.
CC     THE PERIOD IS 2**16=65536, WHICH IS SUFFICIENT FOR THIS
CC     APPLICATION.
       AL=NCAS
       JLP=JLP+1
       IF (JLP.LE.MAXRP) GOTO 10
       RETURN
 10    DO 20 JNC=1,NVAR
 30    JEY=JEY*5761+999
       JEYQT=JEY/65536
       JEY=JEY-JEYQT*65536
       RY=JEY
       RNULE=RY/65536.0
       JRAN=INT(RNULE*AL)+1
       IF (JNC.EQ.1) GOTO 50
       JNCM=JNC-1
       DO 40 J=1,JNCM
       IF (JNDVC(J).EQ.JRAN) GOTO 30
 40    CONTINUE
 50    JNDVC(JNC)=JRAN
 20    CONTINUE
       RETURN
       END
CC
CC  *GENPN* : GENERATES ALL COMBINATIONS OF 'NVAR' INTEGERS
CC            BETWEEN 1 AND 'NCAS', WHERE NVAR IS AT MOST FIVE.
CC
       SUBROUTINE GENPN(NCAS,NVAR,JNDVC,NDXX,JJJ)

C       INCLUDE 'MID_REL_INCL:implicit.inc'
C       DIMENSION JNDVC(NDXX)
       INTEGER J,JJJ,NDXX,NCAS,NVAR
       INTEGER JNDVC(NDXX)

       IF (JJJ.LE.1) THEN
         DO 5 J=1,NVAR
 5         JNDVC(J)=J
           RETURN
       ENDIF
       IF (JNDVC(NVAR).EQ.NCAS) GOTO 10
       JNDVC(NVAR)=JNDVC(NVAR)+1
       RETURN
 10    IF (JNDVC(NVAR-1).EQ.NCAS-1) GOTO 20
       JNDVC(NVAR-1)=JNDVC(NVAR-1)+1
       JNDVC(NVAR)=JNDVC(NVAR-1)+1
       RETURN
 20    IF(JNDVC(NVAR-2).EQ.NCAS-2) GOTO 30
       JNDVC(NVAR-2)=JNDVC(NVAR-2)+1
       JNDVC(NVAR-1)=JNDVC(NVAR-2)+1
       JNDVC(NVAR)=JNDVC(NVAR-1)+1
       RETURN
 30    IF(JNDVC(NVAR-3).EQ.NCAS-3) GOTO 40
       JNDVC(NVAR-3)=JNDVC(NVAR-3)+1
       JNDVC(NVAR-2)=JNDVC(NVAR-3)+1
       JNDVC(NVAR-1)=JNDVC(NVAR-2)+1
       JNDVC(NVAR)=JNDVC(NVAR-1)+1
       RETURN
 40    IF(JNDVC(NVAR-4).EQ.NCAS-4) GOTO 50
       JNDVC(NVAR-4)=JNDVC(NVAR-4)+1
       JNDVC(NVAR-3)=JNDVC(NVAR-4)+1
       JNDVC(NVAR-2)=JNDVC(NVAR-3)+1
       JNDVC(NVAR-1)=JNDVC(NVAR-2)+1
       JNDVC(NVAR)=JNDVC(NVAR-1)+1
       RETURN
 50    IF(JNDVC(NVAR-5).EQ.NCAS-5) RETURN
       JNDVC(NVAR-5)=JNDVC(NVAR-5)+1
       JNDVC(NVAR-4)=JNDVC(NVAR-5)+1
       JNDVC(NVAR-3)=JNDVC(NVAR-4)+1
       JNDVC(NVAR-2)=JNDVC(NVAR-3)+1
       JNDVC(NVAR-1)=JNDVC(NVAR-2)+1
       JNDVC(NVAR)=JNDVC(NVAR-1)+1
       RETURN
       END
CC
CC  *EQUAT* : SOLVES A SYSTEM OF LINEAR EQUATIONS.
CC
      SUBROUTINE EQUAT(AM,NMXV,NDXX,HVEC,JDMB,NA,NB,NERR,DTEM)

C       INCLUDE 'MID_REL_INCL:implicit.inc'
C      DIMENSION HVEC(JDMB),AM(NMXV,NDXX)
C      DOUBLE PRECISION HVEC,TURN,SWAP,DTEM,DETER
      INTEGER J,JBEGC,JBEGX,JDEL,JENDC,JENDX,JDM,JDMB,JHFD
      INTEGER JMAT,JNC,JNCB,JNCC,JNCD,JNCE,JNCF,JNK,JROW
      INTEGER LCLPL,LDEL,N,NA,NB,NC,NDXX,NEQA,NERR,NMXV,NZNDE

      REAL AM(NMXV,NDXX)
      DOUBLE PRECISION HVEC(JDMB),TURN,SWAP,DTEM,DETER

      JDM=NMXV
      DETER=1.0
      N=NA
      JMAT=N+NB
      JNK=0
      DO 10 J=1,JMAT
      JNK=(J-1)*NMXV
      DO 10 NC=1,NMXV
      JNK=JNK+1
      HVEC(JNK)=AM(NC,J)
 10   CONTINUE
      NZNDE=N-1
      LCLPL=-JDM
      DO 120 JHFD=1,N
      TURN=0.
      LCLPL=LCLPL+JDM+1
      JDEL=LCLPL+N-JHFD
      DO 40 JNCB=LCLPL,JDEL
      IF(DABS(HVEC(JNCB))-DABS(TURN)) 40,40,30
 30   TURN=HVEC(JNCB)
      LDEL=JNCB
 40   CONTINUE
      IF(TURN) 50,170,50
 50   IF(LDEL-LCLPL) 60,80,60
 60   DETER=-DETER
      LDEL=LDEL-JDM
      JNCB=LCLPL-JDM
      DO 70 JNCC=JHFD,JMAT
      LDEL=LDEL+JDM
      JNCB=JNCB+JDM
      SWAP=HVEC(JNCB)
      HVEC(JNCB)=HVEC(LDEL)
 70   HVEC(LDEL)=SWAP
 80   DETER=DETER*TURN
      IF (JHFD.EQ.N)  GOTO 120
      TURN=1./TURN
      JNCB=LCLPL+1
      DO 90 JNCC=JNCB,JDEL
 90   HVEC(JNCC)=HVEC(JNCC)*TURN
      JNCD=LCLPL
      JROW=JHFD+1
      DO 110 JNCB=JROW,N
      JNCD=JNCD+1
      JNCE=LCLPL
      JNCF=JNCD
      DO 100 JNCC=JROW,JMAT
      JNCE=JNCE+JDM
      JNCF=JNCF+JDM
 100  HVEC(JNCF)=HVEC(JNCF)-HVEC(JNCE)*HVEC(JNCD)
 110  CONTINUE
 120  CONTINUE
      DTEM=DETER
      NERR=0
      NEQA=N+1
      JBEGX=NZNDE*JDM+1
      DO 150 JNC=NEQA,JMAT
      JBEGX=JBEGX+JDM
      JENDX=JBEGX+N
      JBEGC=N*JDM+1
      JENDC=JBEGC+NZNDE
      DO 140 JNCB=1,NZNDE
      JENDX=JENDX-1
      JBEGC=JBEGC-JDM
      JENDC=JENDC-JDM-1
      HVEC(JENDX)=HVEC(JENDX)/HVEC(JENDC+1)
      SWAP=HVEC(JENDX)
      JNCD=JBEGX-1
      DO 130 JNCC=JBEGC,JENDC
      JNCD=JNCD+1
 130  HVEC(JNCD)=HVEC(JNCD)-HVEC(JNCC)*SWAP
 140  CONTINUE
      HVEC(JBEGX)=HVEC(JBEGX)/HVEC(1)
 150  CONTINUE
      JNC=-JDM
      JBEGX=NZNDE*JDM+1
      JENDX=JBEGX+NZNDE
      DO 160 JNCB=NEQA,JMAT
      JBEGX=JBEGX+JDM
      JENDX=JENDX+JDM
      JNC=JNC+JDM
      JNCD=JNC
      DO 160 JNCC=JBEGX,JENDX
      JNCD=JNCD+1
 160  HVEC(JNCD)=HVEC(JNCC)
      GOTO 180
 170  NERR=-1
      DTEM=DETER
 180  JNK=0
      DO 190 J=1,JMAT
      DO 190 NC=1,NMXV
      JNK=JNK+1
      AM(NC,J)=HVEC(JNK)
 190  CONTINUE
      RETURN
      END
