a34.for

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

      DIMENSION A(ISTOH),B(ISTOH),XM(*),TT(*),W(*),WW(*),R(NN,*),D(*)
      INTEGER MAXA(*),NCOLBV(*),ICOPL(*),NLOC(*)
C
      REWIND NSHIFT
      NJ=1
      BETOL=0.0001D0
      RQTOL=0.0000001D0
      GSTOL=0.001D0
      IF (IFPR.NE.0) WRITE (IOUT,2000)
C
C
   10 IF (NCEV.EQ.0) GO TO 60
      REWIND NT
      READ (NT)
      DO 40 I=1,NCEV
      READ (NT) (TT(K),K=1,NN),(WW(K),K=1,NN)
      AL=0.0D0
      DO 20 K=1,NN
   20 AL=AL + R(K,M1)*WW(K)
      DO 30 K=1,NN
   30 R(K,M1)=R(K,M1) - AL*TT(K)
   40 CONTINUE
C
   60 IF (IMASS.EQ.1) GO TO 65
      CALL MLTPLY (TT,B,R(1,M1),MAXA,NN,NCOLBV,ISTOH,NBLOCK,NMASS)
      GO TO 75
   65 DO 70 K=1,NN
   70 TT(K)=XM(K)*R(K,M1)
C
C
   75 IF (M1.EQ.1) GO TO 110
      IF (KKK.EQ.1) REWIND NSHIFT
      MM=M1 - 1
      DO 90 I=1,MM
      AL=0.0D0
      DO 80 K=1,NN
   80 AL=AL + R(K,M1)*R(K,I)
      DO 85 K=1,NN
   85 TT(K)=TT(K) - AL*R(K,I)
      IF (KKK.GT.1) GO TO 90
      READ (NSHIFT) (WW(K),K=1,NN)
      DO 88 K=1,NN
   88 R(K,M1)=R(K,M1) - AL*WW(K)
   90 CONTINUE
      IF (KKK.EQ.1) GO TO 110
C
      DO 92 K=1,NN
   92 R(K,M1)=TT(K)
      CALL BANDET (A,B,XM,TT,W,MAXA,NCOLBV,ICOPL,NN,ISTOH,NBLOCK,
     1             SHIFT,NSCH,IMASS,FDETA,IDETA,2)
      ALFA=0.0D0
      DO 100 K=1,NN
  100 ALFA=ALFA + TT(K)*R(K,M1)
C
      IF (IMASS.EQ.1) GO TO 105
      CALL MLTPLY (R(1,M1),B,TT,MAXA,NN,NCOLBV,ISTOH,NBLOCK,NMASS)
      GO TO 130
  105 DO 108 K=1,NN
  108 R(K,M1)=XM(K)*TT(K)
      GO TO 130
C
  110 DO 120 K=1,NN
      DUM=TT(K)
      TT(K)=R(K,M1)
  120 R(K,M1)=DUM
      ALFA=0.0D0
C
  130 BETA=0.0D0
      DO 135 K=1,NN
  135 BETA=BETA + TT(K)*R(K,M1)
      D(NJ)=ALFA/BETA
      BETA=SQRT(BETA)
      GAMA=BETA
      DO 140 K=1,NN
      TT(K)=TT(K)/BETA
  140 R(K,M1)=R(K,M1)/BETA
      IF (IFPR.NE.0) WRITE (IOUT,2020) M1,ALFA,BETA,GAMA,D(NJ)
C
      WRITE (NSHIFT) (TT(K),K=1,NN)
      IF (M1.EQ.M2) GO TO 500
      BETA=0.0D0
      DO 210 K=1,NN
  210 WW(K)=0.0D0
      MM=M1 + 1
C
      DO 400 J=MM,M2
      J1=J - 1
C
C     INVERSE ITERATION
C
      DO 220 K=1,NN
  220 R(K,J)=R(K,J1)
      CALL BANDET (A,B,XM,R(1,J),W,MAXA,NCOLBV,ICOPL,NN,ISTOH,NBLOCK,
     1             SHIFT,NSCH,IMASS,FDETA,IDETA,2)
C
C
      ALFA=0.0D0
      DO 230 K=1,NN
  230 ALFA=ALFA + R(K,J)*R(K,J1)
      RQT=ALFA
      DO 240 K=1,NN
  240 R(K,J)=R(K,J) - ALFA*TT(K) - BETA*WW(K)
      DO 250 K=1,NN
  250 WW(K)=TT(K)
C
C
      RQB=ALFA*ALFA + BETA*BETA
      IF (IMASS.EQ.1) GO TO 265
      CALL MLTPLY (TT,B,R(1,J),MAXA,NN,NCOLBV,ISTOH,NBLOCK,NMASS)
      GO TO 275
  265 DO 270 K=1,NN
  270 TT(K)=XM(K)*R(K,J)
  275 BETA=0.0D0
      DO 280 K=1,NN
  280 BETA=BETA + R(K,J)*TT(K)
      RQB=RQB + BETA
      RQ=RQT/RQB
      BETA=SQRT(BETA)
      GAMA=SQRT(RQB)
      DO 285 K=1,NN
  285 R(K,J)=R(K,J)/BETA
      IF (IFPR.NE.0) WRITE (IOUT,2020) J,ALFA,BETA,GAMA,RQ
C
C
      DO 290 I=1,NJ
      IF (ABS(D(NJ) - RQ).LE.ABS(RQ)*RQTOL) GO TO 295
  290 CONTINUE
      IF (BETA/GAMA.GT.BETOL) GO TO 300
  295 DO 298 K=1,NN
  298 R(K,J)=0.0D0
      IJ=NLOC(NCEV + J)
      R(IJ,J)=1.0D0
      M1=J
      NJ=NJ + 1
      GO TO 10
C
C
  300 REWIND NSHIFT
      IFLAG=0
      DO 330 I=1,NJ
      READ (NSHIFT) (TT(K),K=1,NN)
      IF (IFLAG.GT.0) GO TO 330
      AL=0.0D0
      DO 310 K=1,NN
  310 AL=AL + R(K,J)*R(K,I)
      IF (ABS(AL).GT.GSTOL) IFLAG=1
      DO 320 K=1,NN
  320 R(K,J)=R(K,J) - AL*TT(K)
  330 CONTINUE
      IF (IFLAG.GT.0) GO TO 295
C
      DO 340 K=1,NN
  340 TT(K)=R(K,J)
      IF (IMASS.EQ.1) GO TO 350
      CALL MLTPLY (R(1,J),B,TT,MAXA,NN,NCOLBV,ISTOH,NBLOCK,NMASS)
      GO TO 370
  350 DO 360 K=1,NN
  360 R(K,J)=XM(K)*TT(K)
C
  370 AL=0.0D0
      DO 380 K=1,NN
  380 AL=AL + TT(K)*R(K,J)
      AL=SQRT(AL)
      DO 390 K=1,NN
      TT(K)=TT(K)/AL
  390 R(K,J)=R(K,J)/AL
C
      NJ=NJ + 1
      D(NJ)=RQ
      WRITE (NSHIFT) (TT(K),K=1,NN)
  400 CONTINUE
C
C
  500 IF (KKK.GT.1 .OR. IINTER.EQ.0) GO TO 600
      REWIND NSHIFT
      DO 530 J=1,M2
      AL=0.0D0
      DO 510 K=1,NN
  510 AL=AL + R(K,NQ)*R(K,J)
      READ (NSHIFT) (TT(K),K=1,NN )
      DO 520 K=1,NN
  520 R(K,NQ)=R(K,NQ) - AL*TT(K)
  530 CONTINUE
C
  600 RETURN
C
 2000 FORMAT (///,49H STARTING VECTORS ARE GENERATED BY LANCZOS METHOD /
     1     /,7H VECTOR,59H     ALFA           BETA           GAMA
     1  RAYL. QUO.     ,/  )
 2020 FORMAT (I7,4E15.5)
 2050 FORMAT (21H GRAM-SCHMIDT FACTORS,(/8E15.5))
C
      END
      SUBROUTINE RAPID (EIGV,D,TT,W,EVC1,EVC2,RTOLV,R,RR,FREQ,WW,XM,
     1                  NLOC,NSIT,NN,NQ)
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /SOL/ NUMNP,NEQ,NWK,NWM,NWC,NUMEST,MIDEST,MAXEST,NSTE,MAL
      COMMON /EIGIF/ COFQ,RBMSH,IESTYP,NFREQ,NMODE,IFPR,IRBM,
     1               IRSPA,IRSPC,NDIR
      COMMON /SCRAP/ SHIFT,RLQ1,NLQ,NITE,NSTEP,JR,ICONV,NCEV,NP,NQA,IFSS
      COMMON /EL/ IND,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1            ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /DIM/ N0,N1,N2,N3,N4,N5,N6,N7,N8,N9,N10,N11,N12,N13,N14,N15
      COMMON /FREQIF/ ISTOH,N1A,N1B,N1C,N1S
      COMMON /ADDB/ NEQL,NEQR,MLA,NBLOCK
      COMMON /ITELMT/ NSMAX,NITEM,NITEMM,NOVM
      COMMON /FAST/ SHIFT1,SHIFT2,IOVER,IRPC,NSTV,IINTER,JROLD,NSCH1,
     1              IACCN,NJUNK,ISVTYP
      COMMON /DIMSSP/ M3,M4,M5,M6,M7,M8,M9
      COMMON /TAPES/ IIN,IOUT
      COMMON /TAPEIG/ SCALE,NSTIF,NT,NMASS,NRED,NSHIFT,NOVER
      COMMON /TOLS/ RTOL,ALPHA,CTOL,ANORM,RCTOL
C
      REAL A
      COMMON A(1)
      INTEGER IA(1)
      EQUIVALENCE (A(1),IA(1))
C
      DIMENSION EIGV(*),D(*),TT(*),W(*),EVC1(*),EVC2(*),RTOLV(*),R(NN,*)
     1         ,RR(NN,*),FREQ(*),WW(*),XM(*),NLOC(*),NSIT(*)
      DATA NFAC/0/, NGS/0/
C
      RMUS=SHIFT
      SHIFTO=SHIFT
      RITOL=RTOL
      IF (IINTER.GT.0) RITOL=RCTOL
      NP=NFREQ - NCEV
      IF (NQ.LT.NP) RITOL=RCTOL
      IVSHF=0
C
      NEIG=0
      DO 580 I=1,NQ
      IF (RTOLV(I).GT.RITOL) GO TO 590
      NEIG=I
      IF (NLOC(NCEV + I).GT.0) GO TO 580
      NLOC(I + NCEV)=NITE - NSIT(I)
  580 CONTINUE
C
  590 IF (NEIG.EQ.0) GO TO 600
      IF (EIGV(NEIG).GE.SHIFT2) GO TO 595
      IF (NCEV + NEIG.LT.NFREQ) GO TO 600
  595 NFREQ=NCEV + NEIG
      ICONV=1
      JR=0
      IF (IFPR.NE.0) WRITE (IOUT,2060) RTOL,NITE,NFAC,NGS
      GO TO 910
C
  600 IF (NITE.LT.NITEMM) GO TO 700
      IF (IEIG.NE.2) WRITE (IOUT,2070)
      IF (IEIG.EQ.2) WRITE (IOUT,2080) NEIG
      ICONV=2
      IF (IEIG.NE.2) IFSS=0
      IF (IEIG.EQ.2 .AND. NEIG.EQ.0) IFSS=0
      NFREQ=NCEV + NEIG
      JR=0
      GO TO 910
C
C
  700 NC=NFREQ - NCEV
      IF (NC.GT.NQ .OR. IINTER.GT.0) NC=NQ
      IF (NC.GT.0) GO TO 705
      WRITE (IOUT,3010)
      STOP
  705 IF (IACCN.EQ.0) GO TO 1200
      IF (RITOL.EQ.RCTOL) GO TO 720
      NEIG=0
      DO 710 I=1,NQ
      IF (RTOLV(I).GT.RCTOL) GO TO 720
      NEIG=I
  710 CONTINUE
  720 IF (NEIG.LT.NQ) GO TO 730
      IVSHF=1
      GO TO 910
C
  730 IF (NITE.LE.NSTEP - 2) GO TO 910
      DO 740 I=1,NC
      TEMP=D(I) - EVC2(I)
      EVC2(I)=0.0D0
      IF (ABS(TEMP).GT.D(I)*1.E-10) EVC2(I)=(EIGV(I) - D(I))/TEMP
  740 CONTINUE
      IF (IFPR.GT.1)
     1 WRITE (IOUT,2300) (EVC2(I),I=1,NC)
      IF (NITE.EQ.NSTEP - 1) GO TO 910
      IF (IFPR.GT.1)
     1 WRITE (IOUT,2350) (EVC1(I),I=1,NC)
C
C
      DO 750 I=1,NC
      TT(I)=EIGV(NQ)*2.**30
      IF (RTOLV(I).GT.1.D-3 .OR. RTOLV(I).LT.RCTOL) GO TO 750
      IF (EVC2(I).GT.0.99D0 .OR. EVC2(I).LT.1.D-5) GO TO 750
      TEMP=EVC2(I) - EVC1(I)
      TEMP=ABS(TEMP/EVC2(I))
      IF (TEMP.GT.CTOL) GO TO 750
      TT(I)=(EIGV(I) - RMUS)/SQRT(EVC2(I)) + RMUS
      IF (TT(I).LT.RMUS) GO TO 750
      IF (RTOLV(I).GT.RCTOL) NLOC(NCEV + I)=-1
  750 CONTINUE
      IF (IFPR.GT.1)
     1 WRITE (IOUT,2400) (TT(I),I=1,NC)
C
      MLQ=0
      TLQ=0.0D0
      DO 760 I=1,NC
      IF (TT(I).GE.EIGV(NQ)*2.**30) GO TO 760
      IF (TT(I).LT.RMUS) GO TO 760
      MLQ=MLQ + 1
      TLQ=TLQ + TT(I)
  760 CONTINUE
C
      IF (MLQ.NE.0)
     1 RLQ1=(NLQ*RLQ1 + TLQ)/(NLQ + MLQ)
      NLQ=MLQ + NLQ
      IF (IFPR.NE.0)
     1 WRITE (IOUT,2500) RLQ1,NLQ
      IF (RLQ1.GT.0.D0) GO TO 762
      IF (NITE.EQ.NITEMM - 1 .AND. NEIG.GT.0) IVSHF=1
      GO TO 910
C
C
  762 J=1
      DO 763 I=1,NC
      IF (EIGV(I).GT.RMUS) GO TO 765
      J=J + 1
  763 CONTINUE
C
  765 IF (J.GT.NC) J=NC
      NR=J
      DO 770 I=J,NC
      IF (RTOLV(I).GT.RCTOL) GO TO 773
      NR=NR + 1
  770 CONTINUE
C
  773 IF (NR.GT.NC) NR=NC
      K=NR
      BAND=EIGV(1) + (RLQ1 - EIGV(1))/3
C
  775 NR=NR - 1
      IF (NR.LE.J) GO TO 778
      SHIFT=0.5*(EIGV(NR) + EIGV(NR - 1))
      IF (SHIFT.LT.1.01*EIGV(NR - 1)) GO TO 775
      IF (SHIFT.GT.0.99*EIGV(NR)) GO TO 775
      IF (SHIFT.GT.BAND) GO TO 775
      J=NR
      GO TO 790
C
C
  778 IF (J.GT.1) GO TO 788
      IF (NCEV.GT.0) GO TO 780
      SHIFT=0.5*EIGV(1)
      IF (SHIFT.LE.RMUS) GO TO 910
      J=1
      GO TO 790
C
C
  780 NR=NCEV + 2
      FREQ(NCEV + 1)=EIGV(1)
      IF (RTOLV(1).GT.RCTOL) NR=NR - 1
  785 NR=NR - 1
      IF (NR.LT.1) GO TO 910
      SHIFT=0.5*FREQ(NR)
      IF (NR.GT.1) SHIFT=SHIFT + 0.5*FREQ(NR - 1)
      IF (SHIFT.LE.RMUS) GO TO 910
      IF (SHIFT.LT.1.01*FREQ(NR - 1)) GO TO 785
      IF (SHIFT.GT.0.99*FREQ(NR)) GO TO 785
      J=NR - NCEV
      GO TO 790
C
C
  788 TEMP=0.5*(EIGV(K) + EIGV(K - 1))
      IF (TEMP.GT.BAND) IVSHF=1
      GO TO 910
C
C
  790 SI=0.0D0
      II=NC
      DO 800 I=K,NC
      IF (EIGV(I).GE.RLQ1) GO TO 800
      IF (RTOLV(I).LE.RITOL .OR. RTOLV(I).GT.0.01D0) GO TO 800
      TOLI=RITOL/RTOLV(I)
      DI=((EIGV(I) - RMUS)/(RLQ1 - RMUS))**2
      ST = LOG10(TOLI) / LOG10(DI)
      DP=((EIGV(I) - SHIFT)/(RLQ1 - SHIFT))**2
      STP = LOG10(TOLI) / LOG10(DP)
      SII=ST - STP
      IF (SII.LE.SI) GO TO 800
      SI=SII
      II=I
  800 CONTINUE
      I=II
      TOLI=RITOL/RTOLV(I)
      DI=((EIGV(I) - RMUS)/(RLQ1 - RMUS))**2
      ST = LOG10(TOLI) / LOG10(DI)
      DP=((EIGV(I) - SHIFT)/(RLQ1 - SHIFT))**2
      STP = LOG10(TOLI) / LOG10(DP)
C
C
      IF (SI.GT.3.D0) GO TO 820
      NOFAC=0
      NOSI=0
      IF (K.LE.1) GO TO 850
      TEMP=0.5*(EIGV(K) + EIGV(K - 1))
      IF (TEMP.LT.BAND) GO TO 850
      IVSHF=1
      GO TO 910
C
  820 MA=NWK/NN + 1
      NOFAC=(MA*MA + 3*MA)/2
      IF (IMASS.EQ.2) NOFAC=NOFAC + MA
      NP=NQ - JR
      NOSI=(2*MA*NP + 2*NP*NP + 3*NP + 2*NP*JR + 2*NQ*NCEV)*SI
      IF (IMASS.EQ.2) NOSI=NOSI + 2*MA*NP*SI
      NOSI=NOSI + 18*NP*NP*NP*SI/NN
      IF (NOFAC.GE.ALPHA*NOSI) GO TO 850
C
C
      IF (IFPR.NE.0) WRITE (IOUT,2600)
      NSTEP=4 + NITE
      RMUS=SHIFT
      NFAC=NFAC + 1
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
C
      IF (IFPR.NE.0) WRITE (IOUT,2650) NSCH
      IF (NSCH.EQ.NSCH1 + NCEV + J - 1) GO TO 855
      WRITE (IOUT,2680)
      IFSS=0
      GO TO 595
C
  850 IF (IFPR.NE.0) WRITE (IOUT,2690)
  855 IF (IFPR.EQ.0) GO TO 910
      J=J + NCEV + NSCH1
      I=I + NCEV + NSCH1
      WRITE (IOUT,2700) J,I,SHIFT
      WRITE (IOUT,2800) TOLI,DI,DP,ST,STP,NOFAC,NOSI
C
  910 SHIFT=RMUS
      IF (IACCN.EQ.0) GO TO 1200
C
C
      IF (NCEV.EQ.0) GO TO 985
      REWIND NT
      READ (NT)
      BAND=ABS(EIGV(NQ) - SHIFT)
      JJ=JR + 1
      NORTHO=0
      DO 980 J=1,NCEV
      IF (NJUNK.GT.0) GO TO 930
      RT=ABS(FREQ(J) - SHIFT)
      IF (RT.LE.BAND) GO TO 930
      READ (NT)
      GO TO 980
C
  930 READ (NT) (TT(I),I=1,NN),(WW(I),I=1,NN)
      NORTHO=NORTHO + 1
      DO 970 K=JJ,NQA
      AL=0.0D0
      DO 940 I=1,NN
  940 AL=AL + TT(I)*RR(I,K)
      DO 950 I=1,NN
  950 RR(I,K)=RR(I,K) - AL*WW(I)
  970 CONTINUE
  980 CONTINUE
      IF (IFPR.NE.0) WRITE (IOUT,2850) NORTHO
      NGS=NGS + NORTHO
C
C
  985 IF (IOVER.NE.1) GO TO 1000
      IF (NITE.GT.NSTEP-4 .AND. NITE.LE.NSTEP-1) GO TO 1000
      REWIND NOVER
      IF (NJUNK.EQ.0) GO TO 987
      DO 986 J=1,NJUNK
  986 READ (NOVER)
C
  987 JJ=JROLD + 1
      DO 995 J=JJ,NC
      IF (NLOC(NCEV + J).EQ.-1) GO TO 988
      READ (NOVER)
      GO TO 995
C
  988 NLOC(NCEV + J)=0
      READ (NOVER) (TT(K),K=1,NN)
      RC=(EIGV(J) - SHIFTO)/(RLQ1 - SHIFTO)
      RA=1./(1. - RC)
      DO 990 K=1,NN
  990 R(K,J)=TT(K) + (R(K,J) - TT(K))*RA
      K=J + NJUNK
      IF (IFPR.NE.0) WRITE (IOUT,2820) K
  995 CONTINUE
C
C
 1000 IF (IVSHF.EQ.0) GO TO 1200
      K=NEIG
      IF (JR.GT.0) K=JR
      IF (IINTER.GT.0) GO TO 1015
C
      NP=NFREQ - NCEV
      NC=MIN0(2*NP,NP + 8)