C     SIMULATION PROGRAM FOR THREE-DIMENSIONAL MECHANICAL SYSTEMS
C     ***********************************************************
C
C    ----VECNET MAINLINE----
C
      INTEGER*4  STMASS,STNVDR,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,NAME
      COMMON T,TLIMIT,H,SUM(3)
      COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC,
     2 KOUNT
      COMMON/TAG/NTYPE(10),MCON(10,2)
      COMMON/EYCEP/EYE(3),CEN(3),UP(3)
      COMMON/BULK/I3S,LN3,OBJ(6,200)
      COMMON/PAINT/I2S,LN,P(4,200)
      COMMON/LIMITS/XMIN,XMAX,YMIN,YMAX
C<<<<< NETWORK SIZE SPECIFICATION FOLLOWS>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      DIMENSION M(7,10),V(10,3),X(10,3),F(10,3),SIZE(10),SPLENG(10)
      DIMENSION VZERO(10,3),XZERO(10,3),RK(10,3,4),FINALK(10,3)
      DIMENSION NAME(10,2),NOD(7),PRIMEK(10,3)
      DIMENSION PNODE(10,3),SMODEL(7,16,2,2),MODNO(10),IVECT(500)
      DIMENSION DIR(3),PERP(3),DIR2(3)
C  GRAPHIICS INIT
      READ(5,125) XMIN,XMAX,YMIN,YMAX
      READ(5,125)(EYE(I),I=1,3)
      READ(5,125) (CEN(I),I=1,3)
      READ(5,125)(UP(I),I=1,3)
      READ(5,125) DIST
      I2S=1
      I3S=1
      KNT10=0
      IND=0
      CALL ROTER
      CALL CLEAR(.FALSE.)
C    --READ IN AXIS DATA--
C
      READ(5,127)((OBJ(I,J),I=1,6),J=1,3)
127   FORMAT(1X,6F10.4)
C  READ IN SYMBOL VALUES
      DO 50 I=1,7
      READ(5,123) MODNO(I)
      IF(MODNO(I).EQ.0)GO TO 50
      IND=MODNO(I)
      DO 49 IL=1,IND
      READ(5,125) ((SMODEL(I,IL,IE,IC),IC=1,2),IE=1,2)
123   FORMAT(1X,I3)
125   FORMAT(1X,4F10.4)
49    CONTINUE
50    CONTINUE
C   READ IN INDICES FOR ELEMENTS
      READ(5,128) IAX,NN,NB,NMAS,NVDR,NDAM,NSPR,NFOR,BSCALE
128   FORMAT(1X,8I3,F10.4)
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C      --CHECK AND DISPLAY AXIS--
      IST=1
      IF(IAX.EQ.0)GO TO 51
      LN=3
      LN3=3
      CALL MAPPER(DIST)
      CALL VECTOR(1,IVECT,1,LAST,5,0,0)
      CALL DRAW(IVECT,1,LAST)
      PAUSE 1
      IST=LAST+1
51    CONTINUE
C
C<<<< CALCULATE INDEXES FOR DIFFERENT ELEMENTS
      NT=NN-1
      STMASS=1
      STNVDR=NMAS+1
      STDAMP=STNVDR+NVDR
      STSPRG=STDAMP+NDAM
      STFORC=STSPRG+NSPR
      ENDMAS=NMAS
      ENDVDR=ENDMAS+NVDR
      ENDAMP=ENDVDR+NDAM
      ENDSPR=ENDAMP+NSPR
      ENDFOR=ENDSPR+NFOR
C
      NROW=NN
      NCOL=NB
      NSPRNG=NSPR
      NUMASS=NMAS
      CALL TREE(M,SIZE,NOD,NAME,NROW,NCOL)
      CALL READER(V,X,SPLENG,NCOL,NSPRNG,STEP)
C
C<<<<< SET COUNTERS AND STEP SIZE
      KOUNT=0
      H=STEP
      IF(STEP.EQ.0.0) H=(TLIMIT-T)*0.01
      IF(H.EQ.0.0)  GO TO 5001
  100 CONTINUE
      IF((NVDR+NFOR).EQ.0) GO TO 201
      CALL DRIVER(V,X,F,NCOL)
 201  CONTINUE
      CALL VELPOS(V,X,M,NROW,NCOL)
      IF(NDAM.EQ.0) GO TO 203
      CALL DAMPER(V,X,F,SIZE,NROW,NCOL)
 203  CONTINUE
      IF(NSPR.EQ.0) GO TO 204
      CALL SPRING(X,F,SIZE,SPLENG,NCOL,NSPRNG)
 204  CONTINUE
      CALL FORCE1(F,M,NROW,NCOL)
C
C  OUTPUT TO SCREEN
C
C  CALC. COORDS OF ALL NODES.
C
      DO 4000 K=1,3
4000  PNODE(NN,K)=0.0
      LN=0
      DO 4100 I=1,NT
      DO 4100 K=1,3
      PNODE(I,K)=PNODE(MCON(I,1),K)+X(I,K)
4100  CONTINUE
C
C  BEGIN LOOP FOR LINE MATRIX FOR GRAPHIC OUTPUT.
      KNT10=KNT10+1
      IF(KNT10.LT.1)GO TO 888
      KNT10=0
      DO 4800 I=1,NB
      IF(NTYPE(I).EQ.0)GO TO 4800
      IF(NTYPE(I).EQ.9)GO TO 4500
      IF(NTYPE(I).NE.7)GO TO 4998
C    ---DIR. COS. FOR FORCE ARROW.
      R1=F(I,1)*F(I,1)+F(I,2)*F(I,2)+F(I,3)*F(I,3)
      IF(R1.EQ.0.0)GO TO 4800
      R1=SQRT(R1)
      DIR(1)=F(I,1)/R1
      DIR(2)=F(I,2)/R1
      DIR(3)=F(I,3)/R1
      GO TO 4999
C   CALC DIRECTION COSINES.
4998  R1=SQRT(X(I,1)*X(I,1)+X(I,2)*X(I,2)+X(I,3)*X(I,3))
      DIR(1)=X(I,1)/R1
      DIR(2)=X(I,2)/R1
      DIR(3)=X(I,3)/R1
4999  MC=MCON(I,2)
      X2=PNODE(MC,1)-EYE(1)
      Y2=PNODE(MC,2)-EYE(2)
      Z2=PNODE(MC,3)-EYE(3)
      R2=SQRT(X2*X2+Y2*Y2+Z2*Z2)
      ALFA2=X2/R2
      BETA2=Y2/R2
      GAMA2=Z2/R2
      FAC1=DIR(1)*ALFA2+DIR(2)*BETA2+DIR(3)*GAMA2
      FAC2=1./SQRT(1.-FAC1*FAC1)
C  DIR. COS. FOR R3
      PERP(1)=FAC2*(DIR(2)*GAMA2-DIR(3)*BETA2)
      PERP(2)=FAC2*(DIR(3)*ALFA2-DIR(1)*GAMA2)
      PERP(3)=FAC2*(DIR(1)*BETA2-DIR(2)*ALFA2)
CC  PUT END-PTS OF LINES IN MATRIX OBJ.
      IF(NTYPE(I).GT.3)GO TO 4300
C  END-PTS OF LINES FOR VETOR SYMBOLS
      DO 4200 J=1,MODNO(NTYPE(I))
      LN=LN+1
      DO 4200 K=1,3
      OBJ(K,LN)=PNODE(MCON(I,1),K)+X(I,K)*SMODEL(NTYPE(I),J
     1,1,1)+PERP(K)*BSCALE*SMODEL(NTYPE(I),J,1,2)
C
      OBJ(K+3,LN)=PNODE(MCON(I,1),K)+X(I,K)*SMODEL(NTYPE(I),
     1J,2,1)+PERP(K)*BSCALE*SMODEL(NTYPE(I),J,2,2)
4200  CONTINUE
      GO TO 4800
C  END-PTS OF LINES FOR POINT SYMBOLS
4300  CONTINUE
      DIR2(1)=BETA2*PERP(3)-GAMA2*PERP(2)
      DIR2(2)=GAMA2*PERP(1)-ALFA2*PERP(3)
      DIR2(3)=ALFA2*PERP(2)-BETA2*PERP(1)
      DO 4400 J=1,MODNO(NTYPE(I))
      LN=LN+1
      DO 4400 K=1,3
      OBJ(K,LN)=PNODE(MCON(I,2),K)+BSCALE*(DIR2(K)*SMODEL(NTYP
     1E(I),J,1,1)+PERP(K)*SMODEL(NTYPE(I),J,1,2))
C
      OBJ(K+3,LN)=PNODE(MCON(I,2),K)+BSCALE*(DIR2(K)*SMODEL(N
     1TYPE(I),J,2,1)+PERP(K)*SMODEL(NTYPE(I),J,2,2))
4400  CONTINUE
      GO TO 4800
C  END-PTS OF A LINE.
4500  LN=LN+1
      DO 4600 K=1,3
      K1=K+3
      OBJ(K,LN)=PNODE(MCON(I,1),K)
      OBJ(K1,LN)=PNODE(MCON(I,2),K)
4600  CONTINUE
4800  CONTINUE
      LN3=LN
C    ----EYE MOVE----
76    CALL MAPPER(DIST)
C  PASS OBJ TO PERSPECTIVE ROUTINE.
      CALL DRAW(IVECT,IST,0)
      CALL VECTOR(1,IVECT,IST,LAST,7,0,0)
      CALL DRAW(IVECT,IST,LAST)
      CALL PULSE(N,N1,IVECT,IST,LAST)
      IF(IND.EQ.1)GO TO 79
      PAUSE 2
79    IND=1
C
C
C     -----INTEGRATE AND RETURN TO GRAPHICS SECTION----
888   CONTINUE
      CALL RUNGE(M,V,X,F,SIZE,SPLENG,VZERO,XZERO,RK,NROW,NCOL,NSPRNG,
     1 NUMASS,FINALK,PRIMEK)
      IF(T.GT.TLIMIT) GO TO 8
      IF(KOUNT.GT.6000)   GO TO 8
      GO TO 100
C
  8   CONTINUE
      STOP 1
 5001 CONTINUE
      STOP 2
      END
C
      SUBROUTINE READER(V,X,SPLENG,NCOL,NSPRNG,STEP)
C
      INTEGER*4  STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR
      COMMON T,TLIMIT,H,SUM(3)
      COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC,
     2 KOUNT
      DIMENSION V(10,3),X(10,3),SPLENG(10)
C
C<<<< READ INITIAL TIME AND INITIAL BRANCH VALUES
      READ (5,2200) T
 2200 FORMAT (F10.4)
      DO 1 I=STMASS,ENDMAS
      READ(5,2201) J,V(I,1),V(I,2),V(I,3),X(I,1),X(I,2),X(I,3)
 2201 FORMAT( I5,5X,6F10.4)
   1  CONTINUE
C
C<<<< READ UNSTRETCHED SPRING LENGTHS (SPLENG)
      DO 2 I=1,NSPR
      READ(5,2204)  J, SPLENG(I)
 2204 FORMAT(7X,I3,10X,F10.4)
C<<< CHECK IF I AGREES WITH SPRING NUMBER ON DATA CARD
      N=STSPRG+I-1
   2  CONTINUE
C
C<<<< READ INTEGRATION PARAMETERS
      READ(5,2203) TLIMIT,STEP,NUMBER
 2203 FORMAT (2F10.4, 5X,I5)
      RETURN
      END
C     <<<<<<<<<<<<<<<<<<<<<<<
      SUBROUTINE TREE(M,SIZE,NOD,NAME,NROW,NCOL)
C     <<<<<<<<<<<<<<<<<<<<<<<
C
C*****TREE GENERATION SUBROUTINE
      INTEGER*4  STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,NAME,NAME1,NAME2
      COMMON T,TLIMIT,H,SUM(3)
      COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC,
     2 KOUNT
      COMMON/TAG/NTYPE(10),MCON(10,2)
      DIMENSION M(7,10),SIZE(10),NOD(7),NAME(10,2)
C*****SET VALUES OF MATRICES AND COUNTERS
      DO 8 I=1,NN
      DO 8 J=1,NB
      M(I,J)=0
   8  CONTINUE
      DO 1  I=1,NN
      NOD(I)=I
   1  CONTINUE
      IT=0
      IC=NN
C
C*****BEGIN TO INPUT BRANCHES
      DO 4 I=1,NB
      READ(5,2001)J,NAME1,NAME2,NSTART,NEND,NTYPE(I),BRMAG,NTREE
 2001 FORMAT(I5,5X,2A4,2X,3I10,F10.4,I10)
      MCON(I,1)=NSTART
      MCON(I,2)=NEND
      NUMERR=1
C
C*****SORT INTO TREE/COTREE
      IF(IT.EQ.NT) GO TO 2
      IF(NOD(NSTART).EQ.NOD(NEND)) GO TO 2
C
C*****INSERT VALUES IN TREE
      NAME(I,1)=NAME1
      NAME(I,2)=NAME2
      SIZE(I)=BRMAG
      M(NSTART,I)=1
      M(NEND,I)=2
      NSS=NOD(NSTART)
      NNN=NOD(NEND)
      DO 3  I2=1,NN
      IF(NOD(I2).EQ.NSS) NOD(I2)=NNN
   3  CONTINUE
      IT=IT+1
      GO TO 4
C
C*****INSERT VALUES IN COTREE
   2  NAME(IC,1)=NAME1
      NAME(IC,2)=NAME2
      SIZE(IC)=BRMAG
      M(NSTART,IC)=1
      M(NEND,IC)=2
      IC=IC+1
   4  CONTINUE
C
C*****PRINT OUT NAMES OF ELEMENTS IN TREE & COTREE
C
C******BEGIN GENERATING CUTSET MATRIX
      DO 11 I=1,NT
      IF(M(I,I).NE.0.0) GO TO 24
C
C*****CHECK COLUMN ENTRIES BELOW DIAGONAL
      I1=I+1
      IF(I1.GT.NT) GO TO 1000
      DO 12 K3=I1,NT
      IF(M(K3,I).GT.0) GO TO 23
  12  CONTINUE
      NUMERR=2
      GO TO 1000
C
C*****ADD NON-ZERO ROW TO ROW(I)
  23  DO 13 K4=1,NB
      M(I,K4)=M(I,K4)+M(K3,K4)
  13  CONTINUE
C
C*****SET ALL COLUMN VALUES TO ZERO EXCEPT DIAGONAL
  24  DO 14 K1=1,NT
      IF(K1.EQ.I) GO TO 14
  21  IF(M(K1,I).EQ.0) GO TO 14
  25  DO 15 J=1,NB
      M(K1,J)=M(K1,J)+M(I,J)
      IF(M(K1,J).GT.2) M(K1,J)=M(K1,J)-3
  15  CONTINUE
      GO TO 21
  14  CONTINUE
  11  CONTINUE
C
C*****REPLACE TWOS BY +1 OR -1
      NUMERR=3
  26  DO 16 I=1,NT
      IF(M(I,I).EQ.0) GO TO 1000
      IF(M(I,I).EQ.1) GO TO 28
      M(I,I)=1
  27  DO 17 K=NN,NB
      IF(M(I,K).EQ.0) GO TO 17
      IF(M(I,K).EQ.2) GO TO22
      M(I,K)=-1
      GO TO 17
  22  M(I,K)=1
  17  CONTINUE
      GO TO 16
  28  DO 18 K=NN,NB
      IF(M(I,K).GT.1) M(I,K)=-1
  18  CONTINUE
  16  CONTINUE
C
C*****PRINT OUT CUTSET MATRIX
      RETURN
C
C*****ERROR MESSAGES ORIGINATE HERE
 1000 CONTINUE
      STOP 3
      END
C
      SUBROUTINE DRIVER(V,X,F,NCOL)
C     <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      INTEGER*4  STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR
      COMMON T,TLIMIT,H,SUM(3)
      COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC,
     2 KOUNT
      DIMENSION V(10,3),X(10,3),F(10,3)
C
C     THIS SUBROUTINE CONTAINS THE TERMINAL EQUATIONS FOR THE SPECIFIED
C     FUNCTIONS.THERE SHOULD BE ONE EQUATION FOR EACH VECTOR COMPONENT O
C     FUNCTIONS.  IF THERE ARE (NVDR) VELOCITY AND POSITION DRIVERS AND
C     FORCE DRIVERS THEN THERE SHOULD BE (6*NVDR+3*NFOR) EQUATIONS FOLLO
      X(2,1)=20.
      X(2,2)=0.
      X(2,3)=15.
      X(3,1)=0.
      X(3,2)=30.
      X(3,3)=20.
      X(4,1)=30.
      X(4,2)=30.
      X(4,3)=0.
      DO 1000 K=1,3
      V(2,K)=0.
      V(3,K)=0.0
      V(4,K)=0.0
 1000 CONTINUE
      RETURN
      END
C
      SUBROUTINE VELPOS(V,X,M,NROW,NCOL)
      INTEGER*4  STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR
      COMMON T,TLIMIT,H,SUM(3)
      COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC,
     2 KOUNT
      DIMENSION V(10,3),X(10,3),M(7,10)
      ITRANS=1
      N1=1
      N2=NT
      DO 100 I=STDAMP,NB
      IROW=I
      CALL SYSUM (V,M,IROW,N1,N2,ITRANS,NROW,NCOL)
      V(I,1)=-SUM(1)
      V(I,2)=-SUM(2)
      V(I,3)=-SUM(3)
      CALL SYSUM (X,M,IROW,N1,N2,ITRANS,NROW,NCOL)
      X(I,1)=-SUM(1)
      X(I,2)=-SUM(2)
      X(I,3)=-SUM(3)
 100  CONTINUE
      RETURN
      END
C
      SUBROUTINE DAMPER (V,X,F,SIZE,NROW,NCOL)
      INTEGER*4  STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR
      COMMON T,TLIMIT,H,SUM(3)
      COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC,
     2 KOUNT
      DIMENSION UR(3),F(10,3),V(10,3),X(10,3),SIZE(10)
      DO 700 I=STDAMP,ENDAMP
C<<< CALCULATE UNIT VECTOR
      RMAG=X(I,1)*X(I,1)+X(I,2)*X(I,2)+X(I,3)*X(I,3)
      R=SQRT(RMAG)
      IF(R.EQ.0.0) GO TO 702
      UR(1)=X(I,1)/R
      UR(2)=X(I,2)/R
      UR(3)=X(I,3)/R
C<<<< CALCULATE DAMPER FORCE
      D=SIZE(I)*(V(I,1)*UR(1)+V(I,2)*UR(2)+V(I,3)*UR(3))
      F(I,1)=D*UR(1)
      F(I,2)=D*UR(2)
      F(I,3)=D*UR(3)
      GO TO 700
  702 CONTINUE
C<<< PROCEDURE FOR ZERO DAMPER LENGTH
      F(I,1)=SIZE(I)*V(I,1)
      F(I,2)=SIZE(I)*V(I,2)
      F(I,3)=SIZE(I)*V(I,3)
  700 CONTINUE
      RETURN
      END
C
      SUBROUTINE SPRING(X,F,SIZE,SPLENG,NCOL,NSPRNG)
      INTEGER*4  STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR
      COMMON T,TLIMIT,H,SUM(3)
      COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC,
     2 KOUNT
      DIMENSION X(10,3),F(10,3),SIZE(10),SPLENG(10),UR(3)
      N=0
      DO 800 I=STSPRG,ENDSPR
      N=N+1
C<<<CALCULATE UNIT VECTORS
      RMAG=X(I,1)*X(I,1)+X(I,2)*X(I,2)+X(I,3)*X(I,3)
      R=SQRT(RMAG)
      IF(R.EQ.0.0) GO TO 802
      UR(1)=X(I,1)/R
      UR(2)=X(I,2)/R
      UR(3)=X(I,3)/R
C<<< CALCULATE SPRING FORCE
      SPRNGK=SIZE(I)*(R-SPLENG(N))
  803 CONTINUE
      F(I,1)=SPRNGK*UR(1)
      F(I,2)=SPRNGK*UR(2)
      F(I,3)=SPRNGK*UR(3)
      GO TO 800
  802 CONTINUE
C<<< PROCEDURE FOR ZERO SPRING LENGTH
      X1=F(I,1)/SIZE(I)
      X2=F(I,2)/SIZE(I)
      X3=F(I,3)/SIZE(I)
      R= SQRT(X1*X1+X2*X2+X3*X3)
      UR(1)=X1/R
      UR(2)=X2/R
      UR(3)=X3/R
      DIFF=ABS(R-SPLENG(N))
       R=0.0
      SPRNGK=SIZE(I)
      GO TO 803
  800 CONTINUE
      RETURN
      END
C
      SUBROUTINE FORCE1(F,M,NROW,NCOL)
C
C<<< THIS SUBROUTINE CALCULATES THE FORCES IN THE VELOCITY DRIVERS
C
      INTEGER*4  STMASS,STCONM,STDAMP,STSPRG,STFORC,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR
      COMMON T,TLIMIT,H,SUM(3)
      COMMON NN,NB,NT,NVDR,NMAS,NDAM,NCON,NSPR,NFOR,ENDVDR,ENDMAS,
     1 ENDCON,ENDAMP,ENDSPR,ENDFOR,STMASS,STCONM,STDAMP,STSPRG,STFORC,
     2 KOUNT
      DIMENSION F(10,3),M(7,10)
      ITRANS=0
      N1=NN
      N2=NB
      DO 100 I=1,ENDVDR
      IROW=I
      CALL SYSUM(F,M,IROW,N1,N2,ITRANS,NROW,NCOL)
      DO 100 K=1,3
      F(I,K)=-SUM(K)
  100 CONTINUE
  101 CONTINUE
      RETURN
      END
