*     ***********************
      SUBROUTINE MATRIX(SKEW)
*     ***********************
*
*     This subrouting borders the normal matrix with the constraint
*     functions, inverts this enlarged matrix and solves the system.
*
      REAL CF(4),C(51,4)
      INTEGER KS(20)
      INCLUDE 'COMMONS.FOR'
      COMMON/IO/EPS,RELAX,RELAXB,RELAXC,SUMOBS,SUMW,SUMWBS,SW,IPONS,
     *          LMAT,LOUT,NUM
      SAVE
*
      DO 10 I=1,MAXS
      DO 10 J=I+1,MAXS
   10 SMM(J,I)=SMM(I,J)
*
*     Border the matrix and enlarge the vector.
*
      IF(MAXCON.GT.0) THEN
         CALL CONPRE(C,CF)
         DO 20 I=1,MAXCON
         IMAXS=I+MAXS
         V(IMAXS)=CF(I)*SW
         DO 20 J=1,MAXS
         SMM(J,IMAXS)=C(J,I)*SW
   20    SMM(IMAXS,J)=SMM(J,IMAXS)
      END IF
      SKEW=DTSLIN(MAXS+MAXCON)
      OOM=SUMW/(NUM-MAXS+MAXCON)*(SCALE*SCALE)
*
*     Calculate the covariance matrix.
*
      DO 30 I=1,MAXS
      DO 30 J=1,MAXS
   30 SMM(I,J)=SMM(I,J)*OOM
*
*     Test on convergence.
*
      DO 40 I=1,MAXS
      X1=SQRT(SMM(I,I))*EPS
      IF(ABS(V(I)).GT.X1) GO TO 50
   40 CONTINUE
      IF(ICYCLE.EQ.NCYCLE) GO TO 50
      ICYCLE=NCYCLE-1
      FORCED=.TRUE.
*
*     Print correlation matrix.
*
   50 IF(LMAT.NE.1.OR.ICYCLE.NE.NCYCLE) GO TO 90
      WRITE(6,1000)
      IA=1
      LIMIT=MAXS
   60 IF((LIMIT-IA+1).GT.20) LIMIT=IA+19
      WRITE(6,1001) (I,I=IA,LIMIT)
      DO 80 I=1,MAXS
      L=0
      X=SMM(I,I)
      DO 70 J=IA,LIMIT
      L=L+1
      Y=SMM(J,J)*X
   70 KS(L)=100*SMM(I,J)/(SQRT(ABS(Y))*SIGN(1.,Y))+0.5
   80 WRITE(6,1002) I,(KS(J),J=1,L)
      IF(LIMIT.GE.MAXS) GO TO 90
      IA=LIMIT+1
      LIMIT=MAXS
      GO TO 60
   90 RETURN
 1000 FORMAT('0CORRELATION MATRIX=')
 1001 FORMAT('0',4X,20I6)
 1002 FORMAT(I5,20I6)
      END
