a34.for

ADINA84有限元编程学习的好例子

      IF (NQ.LT.NC) GO TO 1010
      GO TO 1200
C
 1010 NC=NQ/2
      IF (NQ.GT.16) NC=NQ - 8
      IF (NP.LT.NC + K) K=NP - NC
      IF (NCEV+K+NQ.GT.NN) K=NN - (NCEV + NQ)
 1015 IF (K.EQ.0) GO TO 1200
      NLQ=0
      RLQ1=0.0D0
      NSTEP=NITE + 4
      NITEMM=NITE + NITEM
      IF (NCEV.NE.0) GO TO 1017
      REWIND NT
      READ (NT)
 1017 CONTINUE
      DO 1050 J=1,K
      FREQ(NCEV + J)=EIGV(J)
      DO 1030 I=1,NN
 1030 TT(I)=(EIGV(J) - SHIFT)*R(I,J)
C
      CALL BANDET (A(N2),A(M3),A(M4),TT,W,A(N1),A(N1A),A(N1B),NN,
     1             ISTOH,NBLOCK,SHIFT,NSCH,IMASS,FDETA,IDETA,2)
C
      IF (IMASS.EQ.1) GO TO 1035
      CALL MLTPLY (R(1,J),A(M3),TT,A(N1),NN,A(N1A),ISTOH,NBLOCK,NMASS)
      GO TO 1040
 1035 DO 1038 I=1,NN
 1038 R(I,J)=XM(I)*TT(I)
 1040 AL=0.0D0
      DO 1042 I=1,NN
 1042 AL=AL + TT(I)*R(I,J)
      AL=SQRT(AL)
      DO 1045 I=1,NN
      TT(I)=TT(I)/AL
 1045 R(I,J)=R(I,J)/AL
      WRITE (NT) (TT(I),I=1,NN),(R(I,J),I=1,NN)
 1050 CONTINUE
C
C
      IF (IFPR.NE.0) WRITE (IOUT,2900) NCEV,NP,K
      NCEV=NCEV + K
      IF (IINTER.EQ.0) GO TO 1060
      IF (NCEV.LE.NFREQ) GO TO 1052
      WRITE (IOUT,3000)
      IVSHF=0
      NCEV=NCEV - K
      GO TO 595
C
C
 1052 IF (NJUNK.EQ.0) GO TO 1060
      NR=NCEV + 1
 1055 NR=NR - 1
      IF (NR.LT.2) GO TO 1060
      RMUS=0.5*(FREQ(NR) + FREQ(NR - 1))
      IF (RMUS.LE.SHIFT) GO TO 1060
      IF (RMUS.LT.1.01*FREQ(NR - 1)) GO TO 1055
      IF (RMUS.GT.0.99*FREQ(NR)) GO TO 1055
      NFAC=NFAC + 1
      SHIFT=RMUS
      IF (IFPR.NE.0) WRITE (IOUT,2920) SHIFT
C
      CALL BANDET (A(N2),A(M3),A(M4),TT,W,A(N1),A(N1A),A(N1B),NN,
     1             ISTOH,NBLOCK,SHIFT,NSCH,IMASS,FDETA,IDETA,1)
C
      IF (IFPR.NE.0) WRITE (IOUT,2650) NSCH
      IF (NSCH.EQ.NSCH1 + NR - 1) GO TO 1060
      WRITE (IOUT,2680)
      IVSHF=0
      NCEV=NCEV - K
      IFSS=0
      GO TO 595
C
 1060 NP=NQ - K
      IF (NP.EQ.0) GO TO 1090
      JR=JR - K
      IF (JR.LT.0) JR=0
      DO 1080 N=1,NP
      J=N + K
      EIGV(N)=EIGV(J)
      NSIT(N)=NSIT(J)
      DO 1070 I=1,NN
 1070 R(I,N)=R(I,J)
 1080 CONTINUE
C
C
 1090 NP=NP + 1
      DO 1092 J=NP,NQ
      EIGV(J)=0.0D0
      NSIT(J)=NITE
 1092 CONTINUE
C
      DO 1097 J=NP,NQ
      IF (NCEV + J.GT.NSTV) GO TO 1098
      READ (NSHIFT) (TT(I),I=1,NN)
      NLOC(NCEV + J)=0
      IF (IMASS.EQ.1) GO TO 1093
      CALL MLTPLY (R(1,J),A(M3),TT,A(N1),NN,A(N1A),ISTOH,NBLOCK,NMASS)
      GO TO 1097
 1093 DO 1095 I=1,NN
 1095 R(I,J)=XM(I)*TT(I)
 1097 CONTINUE
      GO TO 1115
C
 1098 K=J
      DO 1105 J=K,NQ
      DO 1100 I=1,NN
 1100 R(I,J)=0.0D0
      IJ=NLOC(NCEV + NJUNK + J)
      R(IJ,J)=1.0D0
      IF (ISVTYP.GT.0) GO TO 1107
 1105 CONTINUE
      GO TO 1109
C
 1107 M1=NJUNK + J
      CALL STARTV (A(N2),A(M3),A(M4),TT,W,WW,RR,A(N1),A(N1A),A(N1B),
     1             NLOC,D,M1,NQA,NN,ISTOH,NBLOCK,2)
C
 1109 DO 1110 J=K,NQ
 1110 NLOC(NCEV + J)=0
C
 1115 REWIND NT
      READ (NT)
      DO 1150 J=1,NCEV
      READ (NT) (TT(I),I=1,NN),(WW(I),I=1,NN)
      DO 1140 K=NP,NQ
      AL=0.0D0
      DO 1120 I=1,NN
 1120 AL=AL + TT(I)*R(I,K)
      DO 1130 I=1,NN
 1130 R(I,K)=R(I,K) - AL*WW(I)
 1140 CONTINUE
 1150 CONTINUE
C
C
 1200 IF (NBLOCK.EQ.1 .AND. (IACCN.EQ.0.OR.NITE.LT.NSTEP-1)) GO TO 1300
      JROLD=JR
      REWIND NOVER
      IF (NJUNK.EQ.0) GO TO 1204
      DO 1202 J=1,NJUNK
 1202 WRITE (NOVER) (RR(K,J),K=1,NN)
 1204 JJ=JR + 1
      DO 1205 J=JJ,NQ
 1205 WRITE (NOVER) (R(K,J),K=1,NN)
      IF (NBLOCK.EQ.1) GO TO 1300
C
C
      REWIND NOVER
C
 1300 RETURN
C
C
 2060 FORMAT (///,30H CONVERGENCE REACHED FOR RTOL  ,E10.4,2X,
     1        14H  AT ITERATION,I4 /,
     2        37H NUMBER OF FACTORIZATIONS PERFORMED =,I5/,
     3        44H NUMBER OF GRAM-SCHMIDT ORTHOGONALIZATIONS =,I6//)
 2070 FORMAT (1H1,51H*** NO CONVERGENCE IN MAXIMUM NUMBER OF ITERATIONS
     1            ,9HPERMITTED/35H WE ACCEPT CURRENT ITERATION VALUES/
     2            42H THE STURM SEQUENCE CHECK IS NOT PERFORMED  )
 2080 FORMAT (///37HIN LINEARIZED BUCKLING ANALYSIS ONLY ,I5,
     1        26HEIGENVALUES HAVE CONVERGED)
 2300 FORMAT (///44H RATE OF CONVERGENCE ESTIMATES,RI(I + 1) ARE,/(6E15.
     1             5/))
 2350 FORMAT (///44H RATE OF CONVERGENCE ESTIMATES,RI(I    ) ARE,/(6E15.
     1             5/))
 2400 FORMAT (///31H ESTIMATES OF LAMBDA(Q + 1) ARE,/,(6E15.5/))
 2500 FORMAT (///36H AVERAGE ESTIMATE OF LAMDA(Q+1) IS =,E20.10,    /
     1            9H BASED ON,I5,10H ESTIMATES   )
 2600 FORMAT (//51H NOFAC IS LESS THAN ALPHA*NOSI AND SHIFT IS APPLIED/)
 2650 FORMAT (//46H STURM SEQUENCE CHECK PERFORMED AT THIS SHIFT. /
     1          60H ESTIMATED NUMBER OF EIGENVALUES TO THE LEFT OF THIS
     2SHIFT =,I5/)
 2680 FORMAT (14H CHECK FAILED. /
     1        81H ESTIMATED NUMBER OF EIGENVALUES IS MORE THAN THE NUMBE
     2R OF COMPUTED EIGENVALUES. /55H REPEAT SOLUTION WITH A LARGER NUMB
     3ER OF TRIAL VECTORS.   /35H ALSO USE A SMALLER VALUE FOR RTOL //)
 2690 FORMAT (//52H SHIFT CALCULATED DOES NOT SATISFY SHIFTING CRITERIA)
 2700 FORMAT (//,31H SHIFT ESTIMATED TO THE LEFT OF,I3,14HTHE EIGENVALUE
     1,10H BASED ON ,I3,19HTHE EIGENVALUE IS =,E15.5/)
 2800 FORMAT (///11X,4HTOLI,14X,1HD,13X,2HDP,14X,1HT,13X,2HTP,4X,5HNOFAC
     1       ,5X,4HNOSI//,5E15.5,2I9)
 2820 FORMAT (42H OVER-RELAXATION IS PERFORMED FOR VECTOR =,I4)
 2850 FORMAT (//53H GRAM-SCHMIDT ORTHOGONALISATION OF CURRENT TRIAL VECT
     1,36HORS IS PERFORMED W.R.T. THE PREVIOUS,I4,13H EIGENVECTORS  )
 2900 FORMAT (//35H EIGENVECTORS ARE TAKEN OUT        , /
     1          45H NUMBER OF EIGENVECTORS ALREADY REMOVED     =,I5,/,
     2          44H NUMBER OF EIGENPAIRS YET TO BE CALCULATED =,I5/,
     3          49H NUMBER OF EIGENVECTORS CURRENTLY BEING REMOVED =,I5)
 2920 FORMAT (//,81H IN ORDER TO MOVE THE POLE OF ATTRACTION TO THE RIGH
     1T, A SHIFT  IS ALSO PERFORMED /,16H SHIFT APPLIED =,E15.5)
C
 3000 FORMAT (1H1,///,14H *** ERROR ***,/,
     1        55H STORAGE OVERFLOW OCCURED.  IN A GIVEN INTERVAL ONLY A
     2       ,47H MAXIMUM OF FIFTY EIGENVALUES CAN BE CALCULATED  )
 3010 FORMAT (1H1,//45H *** STOP ***, ERROR OCCURED IN EIGENSOLUTION,/,
     159H INCREASE NUMBER OF VECTORS USED (NQ) IN SUBSPACE ITERATION  )
C
      END
      SUBROUTINE SCHECK (EIGV,RTOLV,BUP,BLO,BUPC,NEIV,NC,NEI,RTOL,SHIFT)
C .                                                                   .
C .   P R O G R A M                                                   .
C .                                                                   .
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /EL/ IXY,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1            ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /TAPES/ IIN,IOUT
      COMMON /SCRAP/ RMUSS,RLQ1,NLQ,NITE,NSTEP,JR,ICONV,NCEV,NP,NQ,IFSS
C
      DIMENSION EIGV(*),RTOLV(*),BUP(*),BLO(*),BUPC(*),NEIV(*)
C
      FTOL=0.001D0
      IF (IEIG.EQ.2) FTOL=1.D-07
C
      DO 100 I=1,NC
      BUP(I)=EIGV(I)*(1. + FTOL)
  100 BLO(I)=EIGV(I)*(1. - FTOL)
C
      NROOT=NCEV
      II=NC - NCEV - 1
      DO 120 I=1,II
      IF (RTOLV(I).GT.RTOL) GO TO 130
      IF (BUP(I).LE.BLO(I + 1)) NROOT=NCEV + I
  120 CONTINUE
      IF (RTOLV(NC - NCEV).LE.RTOL) NROOT=NC
  130 IF (NROOT.GE.1) GO TO 200
      WRITE (IOUT,1010)
      STOP
C
C
 200  DO 240 I=1,NROOT
 240  NEIV(I)=1
      IF (NROOT.NE.1) GO TO 260
      BUPC(1)=BUP(1)
      LM=1
      L=1
      I=2
      IF (NC.EQ.1) GO TO 300
      GO TO 295
 260  L=1
      I=2
 270  IF (BUP(I-1).LE.BLO(I)) GO TO 280
      NEIV(L)=NEIV(L)+1
      I=I+1
      IF (I.LE.NROOT) GO TO 270
 280  BUPC(L)=BUP(I-1)
      IF (I.GT.NROOT) GO TO 290
      L=L+1
      I=I + 1
      IF (I.LE.NROOT) GO TO 270
      BUPC(L)=BUP(I-1)
 290  LM=L
      IF (NROOT.EQ.NC) GO TO 300
 295  IF (BUP(I-1).LE.BLO(I)) GO TO 300
      IF (I.GT.NCEV .AND. RTOLV(I - NCEV).GT.RTOL) GO TO 300
      BUPC(L)=BUP(I)
      NEIV(L)=NEIV(L)+1
      NROOT=NROOT+1
      IF (NROOT.EQ.NC) GO TO 300
      I=I+1
      GO TO 295
C
C      FIND SHIFT
C
  300 WRITE (IOUT,1020)
      WRITE (IOUT,1005) (BUPC(I),I=1,LM)
      WRITE (IOUT,1030)
      WRITE (IOUT,1006) (NEIV(I),I=1,LM)
      LL=LM-1
      IF (LM.EQ.1) GO TO 310
 330  DO 320 I=1,LL
 320  NEIV(L)=NEIV(L)+NEIV(I)
      L=L-1
      LL=LL-1
      IF (L.NE.1) GO TO 330
 310  WRITE (IOUT,1040)
      WRITE (IOUT,1006) (NEIV(I),I=1,LM)
      L=0
      DO 340 I=1,LM
      L=L+1
      IF (NEIV(I).GE.NROOT) GO TO 350
 340  CONTINUE
 350  SHIFT=BUPC(L)
      NEI=NEIV(L)
C
      RETURN
C
 1005 FORMAT (1H0,6E22.14)
 1006 FORMAT (1H0,6I22)
 1010 FORMAT (37H0***ERROR   SOLUTION STOP IN *SCHECK*, / 12X,
     1        21HNO EIGENVALUES FOUND., / 1X)
 1020 FORMAT (///,37H UPPER BOUNDS ON EIGENVALUE CLUSTERS  )
 1030 FORMAT (34H0NO OF EIGENVALUES IN EACH CLUSTER  )
 1040 FORMAT (42H0NO OF EIGENVALUES LESS THAN UPPER BOUNDS  )
      END
      SUBROUTINE JACOBI (A,B,X,EIGV,D,N,NWA,RTOL,SHIFT,NSMAX,IFPR)
C .                                                                   .
C .   P R O G R A M                                                   .
C .                                                                   .
C .                                                                   .
C . - - OUTPUT - -                                                    .
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /TAPES/ IIN,IOUT
      DIMENSION A(NWA),B(NWA),X(N,N),EIGV(N),D(N)
C
C
      N1=N + 1
      II=1
      DO 10 I=1,N
      IF (B(II).LT.0.D0) GO TO 3
      IF (A(II).GT.0.D0) GO TO 4
      IF (SHIFT.GT.0.D0) GO TO 4
    3 WRITE(IOUT,2020) II,A(II),B(II)
      STOP
    4 D(I)=A(II)/B(II)
      EIGV(I)=D(I)
   10 II=II + N1 - I
      DO 30 I=1,N
      DO 20 J=1,N
   20 X(I,J)=0.0D0
   30 X(I,I)=1.0D0
      IF (N.EQ.1) GO TO 255
C
C
      NSWEEP=0
      NR=N-1
   40 NSWEEP=NSWEEP+1
      IF(IFPR.EQ.2) WRITE(IOUT,2000)NSWEEP
C
C
      EPS=(.01**NSWEEP)**2
      DO 210 J=1,NR
      JP1=J+1
      JM1=J-1
      LJK=JM1*N - JM1*J/2
      JJ=LJK + J
      DO 210 K=JP1,N
      KP1=K+1
      KM1=K-1
      JK=LJK + K
      KK=KM1*N - KM1*K/2 + K
      EPTOLA=(A(JK)/A(JJ))*(A(JK)/A(KK))
      EPTOLA=ABS(EPTOLA)
      EPTOLB=(B(JK)/B(JJ))*(B(JK)/B(KK))
      IF (EPTOLA.LE.EPS .AND. EPTOLB.LE.EPS) GO TO 210
C
C
      AKK=A(KK)*(B(JK) - B(KK)/A(KK)*A(JK))
      AJJ=A(JJ)*(B(JK) - B(JJ)/A(JJ)*A(JK))
      AB =A(JJ)*(B(KK) - A(KK)/A(JJ)*B(JJ))
      IF (AB.EQ.0.D0) GO TO 43
      CHECK=0.25 + (AKK/AB)*(AJJ/AB)
      GO TO 45
   43 AB=0.5*(ABS(AKK) + ABS(AJJ))
      CHECK=0.0D0
      IF (AB.NE.0.D0) CHECK=(AKK/AB)*(AJJ/AB)
   45 IF ( CHECK ) 50,60,60
   50 WRITE (IOUT,2020) JJ,A(JJ),B(JJ),KK,A(KK),B(KK),JK,A(JK),B(JK)
      STOP
   60 SQCH=ABS(AB)*SQRT(CHECK)
      D1=AB/2.+SQCH
      D2=AB/2.-SQCH
      DEN=D1
      IF(ABS(D2).GT.ABS(D1))DEN=D2
      IF(DEN)80,70,80
   70 CA=0.0D0
      CG=-A(JK)/A(KK)
      GO TO 90
   80 CA=AKK/DEN
      CG=-AJJ/DEN
C
C
   90 IF(N-2)100,190,100
  100 IF(JM1-1)130,110,110
  110 DO 120 I=1,JM1
      IM1=I - 1
      IJ=IM1*N - IM1*I/2 + J
      IK=IM1*N - IM1*I/2 + K
      AJ=A(IJ)
      BJ=B(IJ)
      AK=A(IK)
      BK=B(IK)
      A(IJ)=AJ+CG*AK
      B(IJ)=BJ+CG*BK
      A(IK)=AK+CA*AJ
  120 B(IK)=BK+CA*BJ
  130 IF(KP1-N)140,140,160
  140 LJI=JM1*N - JM1*J/2
      LKI=KM1*N - KM1*K/2
      DO 150 I=KP1,N
      JI=LJI + I
      KI=LKI + I
      AJ=A(JI)
      BJ=B(JI)
      AK=A(KI)
      BK=B(KI)
      A(JI)=AJ+CG*AK
      B(JI)=BJ+CG*BK
      A(KI)=AK+CA*AJ
  150 B(KI)=BK+CA*BJ
  160 IF(JP1-KM1)170,170,190
  170 LJI=JM1*N - JM1*J/2
      DO 180 I=JP1,KM1
      JI=LJI + I
      IM1=I - 1
      IK=IM1*N - IM1*I/2 + K
      AJ=A(JI)
      BJ=B(JI)
      AK=A(IK)
      BK=B(IK)
      A(JI)=AJ+CG*AK
      B(JI)=BJ+CG*BK
      A(IK)=AK+CA*AJ
  180 B(IK)=BK+CA*BJ
  190 AK=A(KK)
      BK=B(KK)
      A(KK)=AK+2.*CA*A(JK)+CA*CA*A(JJ)
      B(KK)=BK+2.*CA*B(JK)+CA*CA*B(JJ)
      A(JJ)=A(JJ)+2.*CG*A(JK)+CG*CG*AK
      B(JJ)=B(JJ)+2.*CG*B(JK)+CG*CG*BK
      A(JK)=0.0D0
      B(JK)=0.0D0
C
C
      DO 200 I=1,N
      XJ=X(I,J)
      XK=X(I,K)
      X(I,J)=XJ+CG*XK
  200 X(I,K)=XK+CA*XJ
  210 CONTINUE
C
C
      II=1
      DO 220 I=1,N
      IF (B(II).LT.0.D0) GO TO 212
      IF (A(II).GT.0.D0) GO TO 215
      IF (SHIFT.GT.0.D0) GO TO 215
  212 WRITE(IOUT,2020) II,A(II),B(II)
      STOP
  215 EIGV(I)=A(II)/B(II)
  220 II=II + N1 - I
      IF(IFPR.LT.2)GO TO 230
      WRITE(IOUT,2030)
      WRITE(IOUT,2010) (EIGV(I),I=1,N)
C
C
  230 DO 240 I=1,N
      TOL=RTOL*ABS(D(I))
      DIF=ABS(EIGV(I)-D(I))
      IF(DIF.GT.TOL)GO TO 280
  240 CONTINUE
C
C
      EPS=RTOL**2
      DO 250 J=1,NR
      JM1=J-1
      JP1=J+1
      LJK=JM1*N - JM1*J/2
      JJ=LJK + J
      DO 250 K=JP1,N
      KM1=K-1
      JK=LJK + K
      KK=KM1*N - KM1*K/2 + K
      EPSA=(A(JK)/A(JJ))*(A(JK)/A(KK))
      EPSA=ABS(EPSA)
      EPSB=(B(JK)/B(JJ))*(B(JK)/B(KK))
      IF (EPSA.LE.EPS .AND. EPSB.LE.EPS) GO TO 250
      GO TO 280
  250 CONTINUE
C
C
  255 II=1
      DO 275 I=1,N
      BB=SQRT(B(II))
      DO 270 K=1,N
  270 X(K,I)=X(K,I)/BB
  275 II=II + N1 - I
      IF (IFPR.GT.0) WRITE (IOUT,2040) NSWEEP
      RETURN
C
C
  280 DO 290 I=1,N
  290 D(I)=EIGV(I)
      IF(NSWEEP.LT.NSMAX)GO TO 40
      WRITE (IOUT,2050)
      STOP
C
 2010 FORMAT(1H0,6E20.12)
 2000 FORMAT(27H0SWEEP NUMBER IN *JACOBI* = ,I4)
 2020 FORMAT (38H0*** ERROR  SOLUTION STOP IN JACOBI    /
     1        31H MATRICES NOT POSITIVE DEFINITE /
     2 (4H II=,I4,6HA(II)=,E20.12,6HB(II)=,E20.12))
 2030 FORMAT(36H0CURRENT EIGENVALUES IN *JACOBI* ARE,/)
 2040 FORMAT (//,33H NUMBER OF SWEEPS IN JACOBI ARE =,I4//)
 2050 FORMAT (1H1,12H ** STOP ** /,
     1            39H NO CONVERGENCE IN *JACOBI ITERATIONS*  /,
     2            26H EIGEN SOLUTION ABANDONED )
C*FILE END
      END