a34.for

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

      SUBROUTINE MSUBSP
C
C
C***ADD:DPR***
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /VAR/ NG,MODEX,IUPDT,KSTEP,ITEMAX,IEQREF,ITE,KPRI,
     1             IREF,IEQUIT,IPRI,KPLOTN,KPLOTE
      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 /DIMSSP/ M3,M4,M5,M6,M7,M8,M9
      COMMON /SEQNCE/ IT1,IT2,IT3,IT4,IT5,IT6,IT7,IT8,IT9,IT10,IT11,
     1                IT12,IT13,IT14,IT15,IT16,IT17,IT18,IT19,IT20,
     1                IT21,IT22,IT23,IT24,IT50,IT56,IT57,IT59,IT60
      COMMON /SOL/ NUMNP,NEQ,NWK,NWM,NWC,NUMEST,MIDEST,MAXEST,NSTE,MA
      COMMON /PRCON/ IDATWR,IPRIC,NPB,IDC,IVC,IAC,IPC,IPNODE(3,15)
      COMMON /TAPEIG/ SCALE,NSTIF,NT,NMASS,NRED,NSHIFT,NOVER
      COMMON /DPR/ ITWO
      COMMON /EIGIF/ COFQ,RBMSH,IESTYP,NFREQ,NMODE,IFPR,IRBM,
     1               IRSPA,IRSPC,NDIR
      COMMON /FREQIF/ ISTOH,N1A,N1B,N1C,N1S
      COMMON /ADDB/ NEQL,NEQR,MLA,NBLOCK
      COMMON /RANDI/ N0A,N1D,IELCPL
      COMMON /DISCON/ NDISCE,NIDM
      COMMON /TOLS/ RTOL,ALPHA,CTOL,ANORM,RCTOL
      COMMON /SCRAP/ SHIFT,RLQ1,NLQ,NITE,NSTEP,JR,ICONV,NCEV,NP,NQ,IFSS
      COMMON /FAST/ SHIFT1,SHIFT2,IOVER,IRPC,NSTV,IINTER,JROLD,NSCH1,
     1              IACCN,NJUNK,ISVTYP
      COMMON /ITELMT/ NSMAX,NITEM,NITEMM,NOVM
      COMMON /TAPES/ IIN,IOUT
      COMMON /CJUNK/ CHEAD
      CHARACTER CHEAD*80
      COMMON A(1)
      REAL A
      DIMENSION IA(1)
      EQUIVALENCE (IA(1),A(1))
C
      IOUT=66
      NWM=NWK
      IF (IMASS.EQ.1) NWM=NEQ
      NMDOF=NEQ
C
      NSTIF = IT4
      IF (KLIN.GT.0) NSTIF = IT12
      IF (IEIG.NE.2) GO TO 25
      NSTIF = IT12
   25 CONTINUE
      NT = IT9
      NRED=IT10
      NMASS = IT11
      NSHIFT = IT18
      NOVER = IT19
C
C
      NREAD=0
      IACCN=0
      IFSS =0
      IRBM =0
      COFQ  =0.0D0
      RBMSH =0.0D0
      SHIFT1=0.0D0
      SHIFT2=0.0D0
      TOL1  =1.D-6
      TOL2  =1.D-10
C
      IF (IEIG.EQ.2) GO TO 28
      IF (IINTER.EQ.1) GO TO 26
C
      READ (IIN,1000) NITEM,IRBM,RBMSH,COFQ,RTOL,ISVTYP,NSTV,IFSS,
     1                IACCN
      GO TO 29
C
   26 READ (IIN,1001) NITEM,IRBM,RBMSH,SHIFT1,SHIFT2,RTOL,ISVTYP,
     1                NSTV
      IACCN=1
      IFSS=1
      GO TO 29
C
   28 READ (IIN,1003) NITEM,RTOL,ISVTYP,NSTV
      IACCN=0
      IFSS=1
   29 CONTINUE
C
      IF (NREAD.GT.0) GO TO 30
C
      M3=N2 + ISTOH*ITWO
      M4=M3
      IF (IMASS.EQ.2 .OR. NBLOCK.GT.1) M4=M3 + ISTOH*ITWO
      M5=M4
      IF (IMASS.EQ.1) M5=M4 + NEQ*ITWO
      GO TO 40
C
C
   30 REWIND NSTIF
      REWIND NMASS
      REWIND NREAD
      READ (NREAD) NEQ,NWK,NWM
      NEQ1=NEQ + 1
      NBLOCK=1
      ISTOH=NWK
      N1A=N1 + NEQ1
      N1B=N1A + 1
      N1C=N1B + 1
      N1D=N1C
      N2=N1D + 2*NBLOCK + 1
      M2=N2
C
      IA(N1A)=NEQ
      IA(N1B)=1
      NBLOC1=2*NBLOCK + 1
      IA(N1D)=0
      IA(N1D + 1)=0
      IA(N1D + 2)=0
C
      M3=M2 + NWK*ITWO
      M4=M3
      IF (IMASS.EQ.2) M4=M3 + NWM*ITWO
      M5=M4
      IF (IMASS.EQ.1) M5=M4 + NWM*ITWO
C
      NN=N1A - 1
      READ (NREAD) (IA(I),I=N1,NN)
      NN=M3 - 1
      READ (NREAD) (A(I),I=M2,NN)
      WRITE (NSTIF) (A(I),I=M2,NN)
      NN=M5 - 1
      READ (NREAD) (A(I),I=M3,NN)
      WRITE (NMASS) (A(I),I=M3,NN)
C
   40 IF (IDATWR.LE.1 .AND. IEIG.NE.2)
     1WRITE (IOUT,2000) CHEAD,NEQ,NWK,NWM
      IF (IDATWR.LE.1 .AND. IEIG.EQ.2)
     1WRITE (IOUT,2001) CHEAD,NEQ,NWK
      IF (MODEX.EQ.0) GO TO 45
C
C
      CALL PNORM (A(N1),A(N1A),A(N2),A(M3),A(M4),IRBM,RBMSH,NEQ,ISTOH,
     1            NBLOCK,NSTIF,NMASS,IMASS,ANORM,NFREQ,NMDOF)
C
   45 NP=MIN0(2*NFREQ,NFREQ + 8)
      IF (NQ.LT.NP) IACCN=1
      IF (NQ.GT.NEQ) NQ=NEQ
      IF (IINTER.EQ.1) GO TO 51
      IF (RTOL.LE.0.D0) RTOL=TOL1
      IF (RTOL.GT.TOL1) write(66,3000)
      GO TO 52
   51 IF (RTOL.LE.0.D0) RTOL=TOL2
      IF (RTOL.LE.TOL2) GO TO 52
      RTOL=TOL2
      write(66,3010)
   52 IF (COFQ.EQ.0.0D0) COFQ=1.D8
      ALPHA=1.0D0
      IOVER=1
      IRPC=0
      IF (IACCN.EQ.1) IRPC=1
      CTOL=0.33D0
      NSMAX=24
      IF (NITEM.EQ.0 .AND. IEIG.NE.2) NITEM=24
      IF (NITEM.EQ.0 .AND. IEIG.EQ.2) NITEM=40
      NITEMM=NITEM
      COFQ2=COFQ*COFQ - RBMSH
      IF (SHIFT2.EQ.0.D0) SHIFT2=COFQ
      PT1=SHIFT1
      SHIFT1=SHIFT1*SHIFT1
      IF (SHIFT1.GT.0.D0) SHIFT1=SHIFT1-RBMSH
      PT2=SHIFT2
      SHIFT2=SHIFT2*SHIFT2 - RBMSH
      NCM=NQ
      IF (NQ.LT.NP) NCM=MIN0(NFREQ + NQ/2 + 1,NFREQ + 8)
      IF (IINTER.EQ.1) NCM=NQ + NFREQ
      IF (NCM.GT.NMDOF) NCM=NMDOF
      NC=NQ
      IF (IDATWR.GT.1) GO TO 70
C
      IF (IEIG.EQ.2) GO TO 75
      IF (IINTER.EQ.1) GO TO 72
C
      write(66,2100)
      write(66,2501) NITEMM
      write(66,2502) IRBM,RBMSH
      write(66,2503) COFQ
      write(66,2505) RTOL,ISVTYP,NSTV
      write(66,2506) IFSS,IACCN
      GO TO 78
C
   72 write(66,2200)
      write(66,2501) NITEMM
      write(66,2502) IRBM,RBMSH
      write(66,2504) PT1,PT2
      write(66,2505) RTOL,ISVTYP,NSTV
      GO TO 78
C
   75 write(66,2300)
      write(66,2501) NITEMM
      write(66,2505) RTOL,ISVTYP,NSTV
   78 CONTINUE
C
C
C
   70 NNC=NC*(NC + 1)/2
      M6=M5 + NNC*ITWO
      M7=M6 + NNC*ITWO
      M8=M7 + NC*NC*ITWO
      M9=M8 + NC*ITWO
      M10=M9 + NC*ITWO
      M11=M10 + NEQ*ITWO
      M12=M11 + NEQ*ITWO
      M13=M12 + NCM
      M14=M13 + NC*ITWO
      M15=M14 + NC*ITWO
      M16=M15 + NC*ITWO
      M17=M16 + NEQ*NC*ITWO
      M18=M17 + NCM*ITWO
      M19=M18 + NEQ*ITWO
      M20=M19 + NC
      CALL SIZE (M20)
C
      IF (MODEX.EQ.0) GO TO 599
C
  100 CALL SSPACE (A(N1),A(N1A),A(N1B),A(N2),A(M3),A(M4),A(M5),A(M6),
     1             A(M7),A(M8),A(M9),A(M10),A(M11),A(M12),A(M13),A(M14),
     2             A(M15),A(M16),A(M17),A(M18),A(M10),A(M11),A(M18),
     3             A(M19),NEQ,NCM,ISTOH,NBLOCK)
C
C
      CALL WRFREQ (A(N2),NFREQ,RBMSH)
C
      IF (NREAD.GT.0) MODEX=0
      IF (NREAD.GT.0) GO TO 599
      M3=N2 + NFREQ*ITWO
      M4=M3 + (NEQ + NDISCE)*ITWO
      M5=M4 + NDOF*NUMNP
      NN=M5 - 1
      REWIND IT9
      READ (IT9) (A(I),I=M4,NN)
C
      CALL WRMOD (A(N2),A(M3),A(M4),NUMNP,NDOF,NEQ,NFREQ,NMODE)
C
  599 CONTINUE
C
      RETURN
C
 1000 FORMAT (2I5,2F10.0,10X,F10.0,4I5)
 1001 FORMAT (2I5,4F10.0,3I5)
 1003 FORMAT (I5,F10.0,2I5)
 2000 FORMAT (1H1,10X,35HS U B S P A C E   I T E R A T I O N,//1X,A//
     155H NUMBER OF EQUATIONS . . . . . . . . . . . . . .(NEQ) =,I8/1X,
     254HNUMBER OF ELEMENTS IN STIFFNESS MATRIX . . . . (NWA) =,I8/1X,
     354HNUMBER OF ELEMENTS IN MASS MATRIX . . . . . . .(NWB) =,I8//)
 2001 FORMAT (1H1,10X,35HS U B S P A C E   I T E R A T I O N,//1X,A//
     155H NUMBER OF EQUATIONS . . . . . . . . . . . . . .(NEQ) =,I8/1X,
     254HNUMBER OF ELEMENTS IN STIFFNESS MATRIX . . . . (NWK) =,I8/)
 2100 FORMAT (////,
     1 59H EIGENVALUE SOLUTION CONTROL DATA                          ,/,
     2 59H FOR LOWEST FREQUENCIES                                      )
 2200 FORMAT (////,
     1 59H EIGENVALUE SOLUTION CONTROL DATA                          ,/,
     2 59H FOR FREQUENCIES WITHIN AN INTERVAL                          )
 2300 FORMAT (////,
     1 59H EIGENVALUE SOLUTION CONTROL DATA                          ,/,
     2 59H FOR LINEARIZED BUCKLING ANALYSIS                            )
C
 2501 FORMAT (//,
     A53H MAXIMUM NUMBER OF SUBSPACE                          ,       /,
     B53H    ITERATIONS ALLOWED . . . . . . . . . . (NITEMM) =,I5      )
 2502 FORMAT (/,
     A53H FLAG FOR PRESENCE OF ZERO FREQUENCIES . . ( IRBM ) =,I5     /,
     B53H    EQ.0, NO ZERO FREQUENCIES                        ,       /,
     C53H    EQ.1, ZERO FREQUENCIES ARE PRESENT               ,      //,
     D53H RIGID BODY MODE SHIFT TO BE APPLIED . . . (RBMSH ) =,E14.6   )
 2503 FORMAT (/,
     A53H CUT-OFF CIRCULAR FREQUENCY  . . . . . . . ( COFQ ) =,E14.6   )
 2504 FORMAT (/,
     A53H CIRCULAR FREQUENCY AT LOWER LIMIT                   ,       /,
     B53H    OF INTERVAL  . . . . . . . . . . . . . (SHIFT1) =,E14.6,//,
     C53H CIRCULAR FREQUENCY AT UPPER LIMIT                   ,       /,
     D53H    OF INTERVAL  . . . . . . . . . . . . . (SHIFT2) =,E14.6   )
 2505 FORMAT (/,
     A53H CONVERGENCE TOLERANCE . . . . . . . . . . (SSTOL ) =,E14.6,//,
     B53H FLAG FOR STARTING VECTOR GENERATION . . . (ISVTYP) =,I5,    /,
     C53H    EQ.0, STANDARD SUBSPACE IS USED                  ,       /,
     D53H    EQ.1, LANCZOS METHOD IS USED                     ,      //,
     E53H NUMBER OF USER-SUPPLIED STARTING                    ,       /,
     F53H    ITERATION VECTORS  . . . . . . . . . . ( NSTV ) =,I5      )
 2506 FORMAT (/,
     A53H FLAG FOR STURM SEQUENCE CHECK . . . . . . ( IFSS ) =,I5     /,
     B53H    EQ.0, NO STURM SEQUENCE CHECK                    ,       /,
     C53H    EQ.1, STURM SEQUENCE CHECK IS APPLIED            ,      //,
     D53H FLAG FOR USAGE OF ACCELERATION SCHEME . . (IACCN ) =,I5     /,
     E53H    EQ.0, NO ACCELERATION                            ,       /,
     F53H    EQ.1, ACCELERATION SCHEME IS USED                         )
C
 3000 FORMAT (///,12H *** WARNING,/,
     1            27H SPECIFIED VALUE FOR RTOL =,E15.5,/,
     2            42H RECOMMENDED VALUE FOR RTOL = .LE. 1.E-06  //)
 3010 FORMAT (///,16H *** WARNING ***,/,
     1 61H RTOL FOR INTERMEDIATE FREQUENCY CALCULATION MUST BE LESS THA,
     2 61HN OR EQUAL TO  1.E-10                                        ,
     3 /,24H RTOL IS RESET TO 1.E-10)
      END
      SUBROUTINE SSPACE (MAXA,NCOLBV,ICOPL,A,B,XM,AR,BR,VEC,EIGV,D,
     1                   TT,W,NLOC,RTOLV,EVC1,EVC2,R,FREQ,WW,BUP,BLO,
     2                   BUPC,NSIT,NN,NCM,ISTOH,NBLOCK)
C .                                                                   .
C .   P R O G R A M                                                   .
C .                                                                   .
C .                                                                   .
C .  - - OUTPUT - -                                                   .
      IMPLICIT DOUBLE PRECISION ( A-H,O-Z )
C***END:DPR***
      COMMON /SOL/ NUMNP,NEQ,NWK,NWM,NWC,NUMEST,MIDEST,MAXEST,NSTE,MAL
      COMMON /EL/ IXY,ICOUNT,NPAR(20),NUMEG,NEGL,NEGNL,IMASS,IDAMP,ISTAT
     1           ,NDOF,KLIN,IEIG,IMASSN,IDAMPN
      COMMON /TAPEIG/ SCALE,NSTIF,NT,NMASS,NRED,NSHIFT,NOVER
      COMMON /TAPES/ IIN,IOUT
      COMMON /FAST/ SHIFT1,SHIFT2,IOVER,IRPC,NSTV,IINTER,JROLD,NSCH1,
     1              IACCN,NJUNK,ISVTYP
      COMMON /SCRAP/ SHIFT,RLQ1,NLQ,NITE,NSTEP,JR,ICONV,NCEV,NP,NQ,IFSS
      COMMON /TOLS/ RTOL,ALPHA,CTOL,ANORM,RCTOL
      COMMON /EIGIF/ COFQ,RBMSH,IESTYP,NFREQ,NMODE,IFPR,IRBM,
     1               IRSPA,IRSPC,NDIR
      COMMON /ITELMT/ NSMAX,NITEM,NITEMM,NOVM
      COMMON /RQSHF/ IRQS
      COMMON /RHSV/ NVEC
      COMMON /BUCK/ SHBU
C
      DIMENSION A(ISTOH),B(ISTOH),XM(*),AR(*),BR(*),VEC(*),EIGV(*),D(*),
     1          TT(*),W(*),RTOLV(*),EVC1(*),EVC2(*),R(NN,*),FREQ(*),
     2          WW(*),BUP(*),BLO(*),BUPC(*),TIM(10)
      INTEGER MAXA(*),NCOLBV(*),ICOPL(*),NLOC(*),NSIT(*)
C
      EQUIVALENCE (NROOT,NFREQ)
C
C
      TOLJ=0.000000000001D0
C
C     INITIALIZATION
C
      ICONV=0
      IRQS=0
      NVEC=1
      NCEV=0
      JR=0
      JROLD=0
      RCTOL=0.0000000001D0
      IF (RTOL.LT.RCTOL) RCTOL=RTOL
      SHIFT=SHIFT1
      NC=NQ
      NNC=NC*(NC + 1)/2
      N1=NC + 1
      NC1=NC - 1
      DO 10 I=1,NC
      NSIT(I)=0
      EVC2(I)=0.0D0
   10 D(I)=0.0D0
      DO 15 I=1,10
   15 TIM(I)=0.0D0
C
C
      REWIND NT
      READ (NT)
      REWIND NSTIF
      REWIND NMASS
      IF (IMASS.EQ.2 .AND. NBLOCK.EQ.1) READ (NMASS) B
      CALL ADTIME ( TIM1 )
C
      ND=NN/NCM
      NEQL=1
      NEQR=0
      MLA=0
      DO 40 NJ=1,NBLOCK
      NCOLB=NCOLBV(NJ)
      NEQR=NEQR + NCOLB
      IF (NBLOCK.GT.1) READ (NSTIF) A
C
      IF (IMASS.EQ.2) GO TO 25
      DO 20 I=NEQL,NEQR
      II=MAXA(I) - MLA
      R(I,1)=XM(I)
      W(I)=XM(I)/A(II)
   20 CONTINUE
      GO TO 35
   25 IF (NBLOCK.GT.1) READ (NMASS) B
      DO 30 I=NEQL,NEQR
      II=MAXA(I) - MLA
      W(I)=B(II)/A(II)
      R(I,1)=B(II)
   30 CONTINUE
C
   35 NEQL=NEQL + NCOLB
      MLA=MAXA(NEQL) - 1
   40 CONTINUE
C
C
      IF (NN.GT.1) GO TO 65
      IF (IMASS.EQ.1) B(1)=XM(1)
      IF (B(1).GT.0.D0) GO TO 62
      WRITE (IOUT,2008)
      STOP
   62 EIGV(1)=A(1)/B(1)
      JR=1
      NSCH=1
      A(1)=1./SQRT(B(1))
      WRITE(NT) A(1)
      NEI=NSCH
      GO TO 1150
C
   65 IF (NCM.EQ.1) GO TO 95
      IF (NC.EQ.1) GO TO 69
      DO 68 J=2,NC
      DO 68 I=1,NN
   68 R(I,J)=0.0D0
   69 L=NN - ND
      DO 90 J=2,NCM
      RT=0.0D0
      DO 70 I=1,L
      IF (W(I).LT.RT) GO TO 70
      RT=W(I)
      IJ=I
   70 CONTINUE
      DO 80 I=L,NN
      IF (W(I).LE.RT) GO TO 80
      RT=W(I)
      IJ=I
   80 CONTINUE
      NLOC(J)=IJ
      W(IJ)=0.0D0
      L=L-ND
      IF (J.LT.NC) R(IJ,J)=1.0D0
   90 CONTINUE
C
      IF (IFPR.EQ.0) GO TO 93
      WRITE (IOUT,2009)
      WRITE (IOUT,2001) (NLOC(J),J=2,NCM)
   93 IF (NC.EQ.1) GO TO 95
C
C
C
      PI=3.141592654D0
      XX=0.5D0
      DO 92 K=1,NN
      XX=(PI + XX)**5
      IX=INT(XX)
      XX=XX - FLOAT(IX)
   92 R(K,NC)=XX
C
C
   95 IF (NSTV.LE.0) GO TO 120
      NV=NSTV
      IF (NV.GT.NC) NV=NC
      REWIND NSHIFT
      READ (NSHIFT)
      DO 110 J=1,NV
      READ (NSHIFT) (TT(I),I=1,NN)
      IF (IMASS.EQ.1) GO TO 100
      CALL MLTPLY (R(1,J),B,TT,MAXA,NN,NCOLBV,ISTOH,NBLOCK,NMASS)
      GO TO 110
  100 DO 105 I=1,NN
  105 R(I,J)=XM(I)*TT(I)
  110 CONTINUE
C
  120 CALL ADTIME ( TIM2 )
      TIM(1)=TIM2 - TIM1
C
C
      CALL BANDET(A,B,XM,TT,W,MAXA,NCOLBV,ICOPL,NN,ISTOH,NBLOCK,SHIFT,
     1            NSCH,IMASS,FDETA,IDETA,1)
C
      CALL ADTIME ( TIM3 )
      TIM(2)=TIM3 - TIM2
      NSCH1=NSCH
      IF (IFPR.NE.0 .AND. IINTER.GT.0)
     * WRITE (IOUT,2060) NSCH1
C
C
      IF (ISVTYP.EQ.0) GO TO 130
      IF (NSTV.GE.NC1 .OR. NC.LE.2) GO TO 130
      M1=NSTV + 1
      DO 125 K=1,NN
  125 R(K,M1)=1.0D0
      CALL STARTV (A,B,XM,TT,W,WW,R,MAXA,NCOLBV,ICOPL,NLOC,EVC1,
     1             M1,NC1,NN,ISTOH,NBLOCK,1)
      CALL ADTIME ( TIM4 )
      TIM(1)=TIM(1) + TIM4 - TIM3
      TIM3=TIM4
C
C
  130 IF (NBLOCK.EQ.1) GO TO 140
      REWIND NOVER
      DO 135 J=1,NC
  135 WRITE (NOVER) (R(K,J),K=1,NN)
      REWIND NOVER
  140 DO 150 J=1,NC
  150 NLOC(J)=0
C
C
      NSTEP=4
      NITE=0
      NLQ=0
      RLQ1=0.0D0
      SHBU=0.0D0
      SHIFB=0.0D0
      NSHBU=0