a34.for

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

      NITESH=0
      BOGTOL = - LOG(RTOL)
  200 NITE=NITE + 1
      IF (IFPR.NE.0) WRITE (IOUT,2010) NITE
C
C
C     SHIFTING CONTROL
C
      IF (IEIG.NE.2) GO TO 219
      IF (ICONV.EQ.1) GO TO 219
      IF (NITE.LT.5) GO TO 216
      IF (NITE.GT.5) GO TO 216
      IF (EIGV(1).LT.0.998D0) GO TO 215
      PCR=500.0D0
      write(66,2150) PCR
      STOP
  215 CONTINUE
      RATIO=EIGV(NQ)/EIGV(1)
      IF (RATIO.GT.1.2D0) GO TO 216
      SHIFB=0.9*EIGV(1)
      GO TO 217
  216 CONTINUE
      NITESH=NITESH+1
      IF (NITE.GE.NITEMM) GO TO 219
      IF (NITESH.LE.4) GO TO 219
      IF (NSHBU.GE.3) GO TO 219
      IF (RTOLV(1).GT.1.D-2) GO TO 219
C
C
      IF (RTOLV(NROOT).LT.TOLJ) GO TO 219
      BOGNEW = - LOG(RTOLV(NROOT))
      DENAB=ABS(BOGNEW - BOGOLD)
      IF (DENAB.LT.TOLJ) GO TO 213
      NAIWS=((BOGTOL-BOGNEW)/(BOGNEW-BOGOLD))*(NITESH-4)*1.2
      NIAVE=NITEMM-NITE+1
      IF (NAIWS.GT.0 .AND. NAIWS.LT.NIAVE) GO TO 219
  213 SHIFB=0.90*EIGV(1)
      SHBUOL=SHBU
  217 SHBU=SHBU+SHIFB
      NITESH=1
      NSHBU=NSHBU+1
      NGOBAN=0
  214 NGOBAN=NGOBAN + 1
      IF (NGOBAN.LT.10) GO TO 212
      write(66,2132)
      STOP
  212 CONTINUE
      CALL BANDET (A,B,XM,TT,W,MAXA,NCOLBV,ICOPL,NN,ISTOH,NBLOCK,
     1             SHBU,NSCH,IMASS,FDETA,IDETA,1)
      IF (NSCH.EQ.0) GO TO 218
      IF (IFPR.NE.0) write(66,2290) NSCH
      SHIFB=0.95*SHIFB
      SHBU=SHBUOL + SHIFB
      GO TO 214
  218 CONTINUE
      IF (IFPR.NE.0) write(66,2300) SHBU
  219 CONTINUE
      IF (IEIG.NE.2 .OR. NSHBU.LT.3) GO TO 202
      IF (NITESH.LE.4) GO TO 202
      IF (RTOLV(NROOT).LT.TOLJ) GO TO 202
      BOGNEW = - LOG(RTOLV(NROOT))
      DENAB=ABS(BOGNEW-BOGOLD)
      IF (DENAB.LT.TOLJ) GO TO 202
      NAIWS=((BOGTOL-BOGNEW)/(BOGNEW-BOGOLD))*(NITESH-4)*0.8
      NIAVE=NITEMM - NITE + 1
      IF (NAIWS.LT.NIAVE) GO TO 202
      write(66,2151) NAIWS,NIAVE
      STOP
C
C
  202 CALL ADTIME ( TIM4 )
      IF (IRPC.EQ.0) JR=0
      JJ=JR + 1
      IF (NBLOCK.EQ.1) GO TO 220
C
C
      NVEC=NC - JR
C
      CALL BANDET (A,B,XM,R(1,JJ),W,MAXA,NCOLBV,ICOPL,NN,ISTOH,NBLOCK,
     1             SHIFT,NSCH,IMASS,FDETA,IDETA,2)
C
      NVEC=1
C
  220 DO 255 J=JJ,NC
      IF (NBLOCK.EQ.1) GO TO 225
      READ (NOVER) (TT(K),K=1,NN)
      GO TO 230
  225 DO 228 K=1,NN
  228 TT(K)=R(K,J)
C
      CALL BANDET (A,B,XM,R(1,J),W,MAXA,NCOLBV,ICOPL,NN,ISTOH,NBLOCK,
     1             SHIFT,NSCH,IMASS,FDETA,IDETA,2)
C
  230 IJ=J
      DO 250 I=1,J
      ART=0.0D0
      IF (I - JR) 238,238,242
  238 DO 240 K=1,NN
  240 ART=ART + R(K,J)*R(K,I)
      ART=ART*(EIGV(I) - SHIFT)
      GO TO 248
  242 DO 246 K=1,NN
  246 ART=ART + R(K,I)*TT(K)
  248 AR(IJ)=ART
  250 IJ=IJ + NC - I
  255 CONTINUE
C
C
      JJ=JR + 1
      DO 310 J=JJ,NC
      IF (IMASS - 1) 275,275,272
  272 CALL MLTPLY (TT,B,R(1,J),MAXA,NN,NCOLBV,ISTOH,NBLOCK,NMASS)
      GO TO 278
  275 DO 276 K=1,NN
  276 TT(K)=XM(K)*R(K,J)
C
  278 IJ=J
      DO 290 I=1,J
      BRT=0.0D0
      IF (ICONV.GT.0) GO TO 279
  277 IF (I - J ) 283,279,279
  279 DO 280 K=1,NN
  280 BRT=BRT + R(K,I)*TT(K)
      GO TO 287
  283 DO 285 K=1,NN
  285 BRT=BRT + R(K,J)*R(K,I)
  287 BR(IJ)=BRT
  290 IJ=IJ + NC - I
      IF (ICONV.GT.0) GO TO 310
      DO 300 K=1,NN
  300 R(K,J)=TT(K)
  310 CONTINUE
      CALL ADTIME ( TIM5 )
      TIM(3)=TIM(3) + TIM5 - TIM4
C
C
      IF (IFPR.NE.2) GO TO 430
      KKK=1
  400 WRITE (IOUT,2020)
      II=1
      DO 410 I=1,NC
      ITEMP=II+NC-I
      WRITE(IOUT,2005) (AR(J),J=II,ITEMP)
  410 II=II + N1 - I
      WRITE (IOUT,2030)
      II=1
      DO 420 I=1,NC
      ITEMP=II+NC-I
      WRITE(IOUT,2005) (BR(J),J=II,ITEMP)
  420 II=II + N1 - I
      IF (KKK - 1) 430,430,440
C
  430 CALL JACOBI (AR,BR,VEC,EIGV,TT,NC,NNC,TOLJ,SHIFT,NSMAX,IFPR)
      CALL ADTIME ( TIM6 )
      TIM(4)=TIM(4) + TIM6-TIM5
C
      IF (IFPR.NE.2) GO TO 440
      WRITE (IOUT,2040)
      KKK=2
      GO TO 400
C
C
  440 DO 445 I=1,NC
  445 EIGV(I)=EIGV(I) + SHIFT
      IF (NC.EQ.1) GO TO 465
C
  448 IS=0
      DO 460 I=1,NC1
      IF (EIGV(I + 1).GE.EIGV(I)) GO TO 460
      IS=IS+1
      EIGVT=EIGV(I + 1)
      EIGV(I + 1)=EIGV(I)
      EIGV(I)=EIGVT
      NCI=NC*I
      NCI1=NC*(I - 1)
      DO 450 K=1,NC
      RT=VEC(NCI + K)
      VEC(NCI + K)=VEC(NCI1 + K)
      VEC(NCI1 + K)=RT
  450 CONTINUE
  460 CONTINUE
      IF (IS.GT.0) GO TO 448
  465 IF (IFPR.EQ.0) GO TO 470
      WRITE (IOUT,2035)
      WRITE (IOUT,2006) (EIGV(I),I=1,NC)
C
C
  470 JRN=0
      IF (NITE.EQ.1 .OR. IACCN.EQ.0 .OR. ICONV.GT.0) GO TO 490
      IF (NC.EQ.1 .OR. EIGV(1).LT.SHIFT1) GO TO 490
      DO 480 I=1,NC1
      DUM=ABS(EIGV(I) - D(I))/EIGV(I)
      IF (DUM.GT.RCTOL) GO TO 490
      IF (1.01*EIGV(I).LT.0.99*EIGV(I + 1)) JRN=I
  480 CONTINUE
  490 IF (JRN.LT.JR) JR=JRN
      IF (JRN.LE.JR) GO TO 500
      JJ=JR + 1
      DO 495 I=JJ,JRN
  495 NLOC(NCEV + I)=NITE - NSIT(I)
C
C
  500 JJ=JR + 1
      DO 540 I=1,NN
      DO 510 J=1,NC
  510 TT(J)=R(I,J)
      DO 530 K=JJ,NC
      KK=NC*(K-1)
      RT=0.0D0
      DO 520 L=1,NC
  520 RT=RT + TT(L)*VEC(KK+L)
  530 R(I,K)=RT
  540 CONTINUE
      IF (IFPR.NE.2) GO TO 542
      WRITE (IOUT,2045)
      DO 541 I=1,NC
      K1=I
      K2=NC*(NC - 1) + I
      WRITE (IOUT,2005) (VEC(J),J=K1,K2,NC)
  541 CONTINUE
C
C
  542 IF (ICONV.GT.0) GO TO 558
      JR=JRN
C
  558 CALL ADTIME ( TIM7 )
      TIM(5)=TIM(5) + TIM7-TIM6
      IF (ICONV.GT.0) GO TO 1000
C
C
      DO 560 I=1,NC
      DIF=ABS(EIGV(I)-D(I))
      DIV=EIGV(I)
      IF (IEIG.EQ.2 .AND. NITESH.EQ.1) DIF=ABS(DIF - SHIFB)
  560 RTOLV(I)=DIF/DIV
      IF (IFPR.NE.0) WRITE (IOUT,2050)
      IF (IEIG.NE.2) GO TO 565
      IF (NITESH.NE.4) GO TO 565
      BOGOLD=20
      IF (RTOLV(NROOT).GT.TOLJ) BOGOLD = - LOG(RTOLV(NROOT))
  565 CONTINUE
      IF (IFPR.NE.0) WRITE (IOUT,2250) (RTOLV(I),I=1,NC)
C
C
  570 CALL ADTIME ( TIM8 )
      NJUNK=0
      DO 590 J=1,NC
      IF (EIGV(J).GT.SHIFT1) GO TO 600
      NJUNK=J
  590 CONTINUE
  600 NJ=NJUNK + 1
      NR=NC - NJUNK
C
      CALL RAPID (EIGV(NJ),D(NJ),TT,W,EVC1(NJ),EVC2(NJ),RTOLV(NJ),
     1            R(1,NJ),R,FREQ,WW,XM,NLOC,NSIT(NJ),NN,NR)
      IF (IEIG.EQ.2 .AND. NITESH.LT.4) ICONV=0
C
      CALL ADTIME ( TIM9 )
      TIM(6)=TIM(6) + TIM9-TIM8
C
  910 DO 920 I=1,NC
      EVC1(I)=EVC2(I)
      EVC2(I)=D(I)
  920 D(I)=EIGV(I)
      IF (JR.EQ.0) GO TO 960
      IF (ICONV.GT.0) GO TO 960
C
C
      II=1
      DO 940 I=1,JR
      IJ=II
      AR(IJ)=EIGV(I) - SHIFT
      BR(IJ)=1.0D0
      IJ=IJ + 1
      IF (I.EQ.JR) GO TO 940
      JJ=I + 1
      DO 930 J=JJ,JR
      AR(IJ)=0.0D0
      BR(IJ)=0.0D0
  930 IJ=IJ + 1
  940 II=II + N1 - I
C
  960 GO TO 200
C
C
 1000 CALL ADTIME ( TIM10 )
      JR=NROOT
      J=0
      DO 1010 I=1,NC
      IF (EIGV(I).LT.SHIFT1) GO TO 1010
      J=J + 1
      FREQ(NCEV + J)=EIGV(I)
 1010 CONTINUE
C
      IF (NROOT.EQ.0) GO TO 1160
      REWIND NT
      READ (NT)
      IF (NCEV.EQ.0) GO TO 1030
      DO 1020 L=1,NCEV
 1020 READ (NT)
 1030 NR=NROOT - NCEV
      IF (NR.EQ.0) GO TO 1080
C
      J=0
      DO 1070 L=1,NC
      IF (EIGV(L).LT.SHIFT1) GO TO 1070
      J=J + 1
      IF (J.GT.NR) GO TO 1080
      IF (IMASS - 1) 1050,1050,1040
 1040 CALL MLTPLY ( W,B,R(1,L),MAXA,NN,NCOLBV,ISTOH,NBLOCK,NMASS)
      GO TO 1060
 1050 DO 1055 I=1,NN
 1055 W(I)=XM(I)*R(I,L)
 1060 WRITE (NT) (R(I,L),I=1,NN),(W(I),I=1,NN)
 1070 CONTINUE
 1080 REWIND NT
      READ (NT)
C
C
      REWIND NSTIF
      IF (NBLOCK.EQ.1) READ (NSTIF) A
      DO 1140 L=1,NROOT
      RT=FREQ(L)
      IF (IEIG.EQ.2) RT=RT + SHBU
      READ (NT) (WW(I),I=1,NN),(R(I,1),I=1,NN)
      CALL MLTPLY (TT,A,WW,MAXA,NN,NCOLBV,ISTOH,NBLOCK,NSTIF)
      VNORM=0.0D0
      IF (IEIG.NE.2) GO TO 1088
      DO 1081 I=1,NN
 1081 VNORM=VNORM + R(I,1)*R(I,1)
      WNORM=0.0D0
      BINFU=1. + TOLJ
      BINFL=1. - TOLJ
      IF (RT.GT.BINFU .OR. RT.LT.BINFL) GO TO 1082
      W(L)=1.D10
      GO TO 1140
 1082 CONTINUE
      PCRIT=1./(1.-RT)
      DO 1083 I=1,NN
      R(I,1)=R(I,1) - PCRIT*(R(I,1)-TT(I))
 1083 WNORM=WNORM + R(I,1)*R(I,1)
      VNORM=SQRT(VNORM)
      WNORM=SQRT(WNORM)
      W(L)=WNORM/VNORM
      GO TO 1140
 1088 CONTINUE
      DO 1100 I=1,NN
 1100 VNORM=VNORM + TT(I)*TT(I)
      WNORM=0.0D0
      DO 1120 I=1,NN
      TT(I)=TT(I) - RT*R(I,1)
 1120 WNORM=WNORM + TT(I)*TT(I)
      VNORM=SQRT(VNORM)
      WNORM=SQRT(WNORM)
      W(L)=WNORM/VNORM
 1140 CONTINUE
C
      WRITE (NT) (FREQ(I),I=1,NROOT)
C
      IF (IEIG.NE.2) WRITE (IOUT,2100)
      IF (IEIG.EQ.2) WRITE (IOUT,2101) SHBU
      IF (RBMSH.NE.0.D0) write(66,2105) RBMSH
      WRITE (IOUT,2006) (FREQ(I),I=1,NROOT)
      WRITE (IOUT,2110) (NLOC(I),I=1,NROOT)
      IF (IEIG.NE.2) WRITE (IOUT,2115)
      IF (IEIG.EQ.2) WRITE (IOUT,2116)
      WRITE (IOUT,2006) (W(I),I=1,NROOT)
      IF (IEIG.NE.2) GO TO 1250
      IPRNOT=0
      DO 1210 I=1,NROOT
      IF (W(I).GT.(1.D-05)) IPRNOT=IPRNOT+1
 1210 CONTINUE
      IF (IPRNOT.NE.0) write(66,2400)
 1250 CONTINUE
      IF (IFSS.EQ.0) GO TO 1160
C
C
      NEI=NROOT
      NJUNK=0
      DO 1142 L=1,NC
      IF (EIGV(L).GT.SHIFT1) GO TO 1145
      NJUNK=L
 1142 CONTINUE
 1145 NJ=NJUNK + 1
      NCM=NCEV + NC - NJUNK
      CALL SCHECK (FREQ,RTOLV(NJ),BUP,BLO,BUPC,NLOC,NCM,NEI,RTOL,SHIFT)
      IF (IEIG.EQ.2) SHIFT=SHIFT + SHBU
C
      WRITE (IOUT,2120) SHIFT
C
C
      CALL BANDET(A,B,XM,TT,W,MAXA,NCOLBV,ICOPL,NN,ISTOH,NBLOCK,SHIFT,
     1            NSCH,IMASS,FDETA,IDETA,1)
C
      NSCH=NSCH - NSCH1
 1150 WRITE (IOUT,2130) NSCH,SHIFT1,SHIFT,NEI
      IF (NSCH.EQ.NEI) GO TO 1160
      IF (IEIG.NE.2) write(66,2132)
      IF (IEIG.NE.2 .AND. NSCH.GT.NEI) WRITE (IOUT,2140)
      IF (IEIG.EQ.2) PCR=1./(1.-EIGV(1)-SHBU)
      IF (IEIG.EQ.2 .AND. NSCH.GT.NEI) WRITE (IOUT,2150) PCR
      STOP
 1160 CALL ADTIME ( TIM11 )
      TIM(7)=TIM(7) + TIM11-TIM10
      TIM(10)=TIM11 - TIM1
      IF (IFPR.GT.0)
     *  WRITE (IOUT,2200) NITE,(TIM(I),I=1,7),TIM(10)
C
      RETURN
C
 2001 FORMAT (1H0,10I10)
 2005 FORMAT (1H ,12E11.4)
 2006 FORMAT (1H0,6E22.14)
 2008 FORMAT (38H0***ERROR   SOLUTION STOP IN 'SSPACE'   , / 12X,
     1        23HNO EIGENVALUES COMPUTED, / 1X)
 2009 FORMAT ( ///,62H DEGREES OF FREEDOM EXCITED BY UNIT STARTING ITERA
     1TION VECTORS)
 2010 FORMAT (1H1,32HI T E R A T I O N   N U M B E R ,I4//)
 2020 FORMAT (28H0PROJECTION OF A (MATRIX AR) )
 2035 FORMAT (30H0EIGENVALUES OF AR-LAMBDA*BR  )
 2030 FORMAT (28H0PROJECTION OF B (MATRIX BR) )
 2040 FORMAT (40H0AR AND BR AFTER JACOBI DIAGONALIZATION  )
 2045 FORMAT (29H0Q MATRIX                       /)
 2050 FORMAT (43H0RELATIVE TOLERANCE REACHED ON EIGENVALUES  )
 2060 FORMAT (//,37H NUMBER OF EIGENVALUES BELOW SHIFT1 =,I5)
 2100 FORMAT (1H1,22H E I G E N V A L U E S )
 2101 FORMAT (1H1,73H I N F O R M A T I O N   O N   T H E   G A M M A
     1 E I G E N V A L U E S  ,//,
     223H EIGENVALUES SHIFTED BY,E20.12,14H ARE EQUAL TO  )
 2105 FORMAT (20H (SHIFTED BY RBMSH =,E20.12,2H ) )
 2110 FORMAT (///59H NUMBER OF SUBSPACE ITERATIONS PERFORMED FOR EACH EI
     1GENPAIR/,(1H0,20I5))
 2115 FORMAT (//,1X,45HPRINT PHYSICAL ERROR NORMS ON THE EIGENVALUES  )
 2116 FORMAT (//,1X,98HTHE PHYSICAL ERROR NORMS ON THE LAMBDA EIGENVALUE
     1S CALCULATED FROM THE ABOVE GAMMA EIGENVALUES ARE )
 2120 FORMAT  (///,23H CHECK APPLIED AT SHIFT  ,E22.14)
 2130 FORMAT (/40H BASED ON STURM SEQUENCE CHECK THERE ARE,I4,
     1         29H EIGENVALUES BETWEEN SHIFT1 =,E15.5,12H AND SHIFT =,
     2         E15.5,//46H NUMBER OF EIGENVALUES CALCULATED BY THE PROGR
     3        ,21HAM IN THIS INTERVAL =,I4 //)
 2132 FORMAT (//,44H *** REQUESTED EIGENVALUES NOT OBTAINED ***  ,/,
     1          50H *** STOP OF SOLUTION ***                          //
     2          50H FAILURE OF SOLUTION ALGORITHM CAN BE DUE TO USE   /
     3          47H OF BAD MODEL OR INAPPROPRIATE USE OF SOLUTION
     4          ,15HPARAMETERS      //)
 2140 FORMAT (95H TO CALCULATE THE MISSING EIGENVALUES REPEAT THE SOLUTI
     1ON USING LARGER NUMBER OF TRIAL VECTORS.,/19H ALSO DECREASE RTOL )
 2150 FORMAT (//,5X,55HTHE SOLUTION OF THE LINEARIZED BUCKLING ANALYSIS
     1FAILED,/,
     25X,48HUPPER BOUND FOR FIRST CRITICAL LOAD OF THE MODEL,/,
     25X,E12.5,19H X (LOAD AT STEP 1),/,
     35X,31HFOR OBTAINING THE EIGENSOLUTION,/,
     410X,33HA) EITHER USE A TIGHTER TOLERANCE,/,
     510X,40HB) OR INCREASE THE LOAD LEVEL FOR STEP 1,/,
     610X,20HC) USE 2*NQ IN MCC 5 )
 2151 FORMAT (5X,54HNUMBER OF ITERATIONS STILL REQUIRED FOR CONVERGENCE
     1 =,I5,/,
     2        5X,54HNUMBER OF ITERATIONS AVAILABLE
     3 =,I5,/,
     45X,63HRERUN ALLOWING FOR MORE ITERATIONS AND/OR USING NQ=3*NQ (MCC
     5 5))
 2200 FORMAT (1H1,///23H EIGENSOLUTION TIME LOG ,///1X,
     A51HNUMBER OF SUBSPACE ITERATIONS PERFORMED           =,I5,/1X,
     151HTIME FOR CALCULATION OF STARTING SUBSPACE         =,F9.3,/1X,
     251HTIME FOR LDLT FACTORIZATION OF STIFFNESS MATRIX   =,F9.3,/1X,
     351HTIME FOR CALCULATION OF PROJECTIONS OF A AND B    =,F9.3,/1X,
     451HTIME FOR SOLVING EIGENSYSTEM OF SUBSPACE OPERATORS=,F9.3,/1X,
     551HTIME FOR SORTING EIGENVALUES, NORMALISING VECTORS =,F9.3,/1X,
     651HTIME FOR CALCULATING AND APPLYING SHIFTS          =,F9.3,/1X,
     751HTIME FOR ERROR NORMS AND STURM SEQUENCE CHECK     =,F9.3,//1X,
     851HTIME FOR EIGENSOLUTION                            =,F9.3,//)
 2250 FORMAT (6E15.5/)
 2290 FORMAT (/,5X,25HWHEN SHIFTING JUMPED OVER,I5,11HEIGENVALUES,/,
     1          5X,20HTHE SHIFT IS REDUCED)
 2300 FORMAT (//,5X,49HSHIFT IMPOSED ON EIGENVALUE PROBLEM FOR BUCKLING=
     1,E15.5)
 2400 FORMAT (//,20H * *   WARNING   * * ,/,
     15X,44HPHYSICAL ERROR NORMS ARE GREATER THAN 1.D-05 ,/,
     25X,65HTO OBTAIN A MORE ACCURATE SOLUTION TIGHTEN THE TOLERANCE (SS
     3TOL)   ,//)
C
      END
      SUBROUTINE STARTV (A,B,XM,TT,W,WW,R,MAXA,NCOLBV,ICOPL,NLOC,D,
     1                   M1,M2,NN,ISTOH,NBLOCK,KKK)
C
C .                                                                   .
C .                                                                   .
C
C
      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 /EIGIF/ COFQ,RBMSH,IESTYP,NFREQ,NMODE,IFPR,IRBM,
     1               IRSPA,IRSPC,NDIR
      COMMON /RQSHF/ IRQS
      COMMON /RHSV/ NVEC
C