subspc.for

著名的berkeley大学开发的结构动力非线性分析程序

C*****
      SUBROUTINE SUBSPC
     $ (LL,AKA,AMA,AV,W,N,MT,MW,VV,SKS,SMS,SV,Y,ID,NMAX,EPS)
C**   LL(N) = INPUT VARIABLE BAND MATRIX STORED INDEX
C**   AKA(NN) = INPUT STIFFNESS MATRIX
C**   AMA(NN) = INPUT MASS MATRIX (IF ID=1)
C**   AMA(N) = INPUT DIAGONAL MASS (IF ID=0)
C**   AV(N,MT) = OUTPUT EIGEN VECTORS
C**   W(MT) = EIGENVALUES (OMAGA SQUARES)
C**   VV(N) = WORKING ARRAY
C**   SKS(MT*(MT+1)/2) = WORKING ARRY
C**   SMS(MT*(MT+1)/2) = WORKING ARRY
C**   SV(MT,MT) = WORKING ARRAY
C**   Y(MT) = WORKING ARRAY
C**   N = ORDER OF AKA AND AMA
C**   MT = TOTAL NO. OF EIGENVECTORS USED IN ITERATON
C**   MW = REQD. NO. OF EIGENVECTORS
C**   ID = 0, AMA IS A DIG MATRIX
C**      = 1, AMA IS STORED AS AKA
C**   NMAX = ALLOWED MAX. NO. OF SUBSPACE ITERATIONS
C**   EPS = TOLERENCES OF EIGENVALES W(MW)
C**   NN = LL(N)+N-1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION LL(1),AKA(1),AMA(1),AV(1,1),W(1)
      DIMENSION VV(1),SKS(1),SMS(1),SV(1),Y(1)
C
      DO 1 I=1,N
      IF (AMA(I).GT.0.) THEN
          AMAMIN=AMA(I)
          GOTO 2
      ENDIF
    1 CONTINUE
    2 DO 5 I=1,N
      IF(AMA(I).LE.0.) GOTO 5
      AMAMIN=DMIN1(DABS(AMA(I)),AMAMIN)
    5 LL(I)=LL(I)-I+1
      AMAMIN=AMAMIN/10000.
      DO 8 I=1,N
    8 IF (AMA(I).LE.0.0) AMA(I)=DABS(AMAMIN)
      NPT=1
      M=MT
      MMM=NMAX
      NNN=0
      KKK=0
      CALL VBDECP(AKA,LL,N)
      CALL SETMVO(LL,AKA,AMA,AV,N,M,ID,VV)
      DO 20 J=1,M
   20 W(J)=0.0
   30 NNN=NNN+1
      DO 50 J=1,M
      Y(J)=W(J)
      JV=(J-1)*N+1
      KK=(J-1)*J/2+1
      DO 40 I=1,N
   40 VV(I)=AV(I,JV)
      CALL VBSOLX(AKA,AV(1,JV),LL,N)
   50 CALL CMPYAT(AV,VV,SKS(KK),J,N)
      IF(NNN/NPT*NPT.NE.NNN.AND.NNN.NE.MMM) GO TO 60
C**   AV  =AFA*AMA*V0=V1
C**   SKS=VIT*AKA*V1
   60 DO 80 J=1,M
      JV=(J-1)*N+1
      MM=(J-1)*J/2+1
      IF(NNN-MMM) 65,75,75
   65 DO 70 I=1,N
   70 VV(I)=AV(I,JV)
      CALL VBCMPY(AMA,VV,AV(1,JV),LL,N,ID)
      GO TO 80
   75 CALL VBCMPY(AMA,AV(1,JV),VV,LL,N,ID)
   80 CALL CMPYAT(AV,VV,SMS(MM),J,N)
      IF(NNN/NPT*NPT.NE.NNN.AND.NNN.NE.MMM) GO TO 90
C**   AV=AMA*V1.OR.AV=V1 (NNN.GE.MMM)
C**   SMS=VIT*AMA*V1
   90 CALL BVEXAV(SMS,SKS,SV,M,W,VV,KKK,1.0D-5,M)
      CALL ORDERX(W,SV,VV,M)
      CALL AMPY  (AV,SV,VV,N,M)
C**   AV  =AMA*V1*SV=AMA*V2.OR.AV=V2 (NNN.GE.MMM)
      IF(NNN-MMM) 110,200,300
  110 DO 120 J=1,MW
      IF(ABS(Y(J)-W(J))-EPS*W(J)) 120,120,30
  120 CONTINUE
      MMM=0
      GO TO 30
  200 CONTINUE
  300 CONTINUE
      RETURN
      END
C*****
      SUBROUTINE SETMVO(LL,AKA,AMA,AV,N,M,ID,VV)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION LL(1),AKA(1),AMA(1),AV(1,1),VV(1)
      DATA IX/1/
      DO 15 J=1,M
      JV=(J-1)*N+1
      DO 10 I=1,N
   10 AV(I,JV)=RANDOM(IX)
   15 AV(J,JV)=AV(J,JV)+1.0
      RETURN
      END
C*****
      FUNCTION RANDOM(IX)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C**   RANDOM NUMBER BY LINEAR CONGRUENTIAL METHOD
C**   FOR 32 BITS INTEGER*4 51603,2147483647
C**   FOR 16 BITS INTEGER*2 403,32767
C**   FOR 1'S COMPLEMENT, CANCEL +1
      IX=51603*IX+3
      IF(IX.LT.0) IX=IX+2147483647+1
   20 RANDOM=IX/2147483647.0
      RETURN
      END
C*****
      SUBROUTINE VBCMPY(A,X,Y,LL,N,ID)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(1,1),X(1),Y(1),LL(1)
      JA=1
      IF(ID.EQ.0) GO TO 40
      DO 20 J=1,N
      YJ=0.0
      XJ=X(J)
      J1=J-1
      JM=J-LL(J)+JA
      JA=LL(J)
      IF(JM.GT.J1) GO TO 20
      DO 10 I=JM,J1
      Y(I)=Y(I)+A(I,JA)*XJ
   10 YJ=YJ+A(I,JA)*X(I)
   20 Y(J)=YJ+A(J,JA)*XJ
      RETURN
   40 DO 50 I=1,N
   50 Y(I)=A(I,JA)*X(I)
      RETURN
      END
C*****
      SUBROUTINE CMPYAT(A,B,C,M,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(1,1),B(1),C(1)
      DO 20 J=1,M
      JA=(J-1)*N+1
      CJ=0.0
      DO 10 I=1,N
   10 CJ=CJ+A(I,JA)*B(I)
   20 C(J)=CJ
      RETURN
      END
C*****
      SUBROUTINE AMPY(A,B,AI,N,M)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(1,1),B(1,1),AI(1)
      DO 30 I=1,N
      KA=1
      DO 10 K=1,M
      AI(K)=A(I,KA)
   10 KA=KA+N
      JA=1
      JB=1
      DO 30 J=1,M
      AIJ=0.0
      DO 20 K=1,M
   20 AIJ=AIJ+AI(K)*B(K,JB)
      A(I,JA)=AIJ
      JA=JA+N
   30 JB=JB+M
      RETURN
      END
C*****
      SUBROUTINE ORDERX(X,V,L,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(1),V(1,1),L(1)
      DO 10 I=1,N
   10 L(I)=I
      M=N
   20 NE=M-1
      IF(NE) 60,60,30
   30 DO 50 I=1,NE
      J=L(I)
      K=L(I+1)
      IF(X(J)-X(K)) 50,50,40
   40 L(I)=K
      L(I+1)=J
      M=I
   50 CONTINUE
      IF(M-NE) 20,20,60
   60 N1=N-1
      DO 90 J=1,N1
      I=J
   70 I=L(I)
      IF(I-J) 70,90,80
   80 II=(I-1)*N+1
      JJ=(J-1)*N+1
      T=X(I)
      X(I)=X(J)
      X(J)=T
      DO 85 KK=1,N
      T=V(KK,II)
      V(KK,II)=V(KK,JJ)
   85 V(KK,JJ)=T
   90 CONTINUE
      RETURN
      END
C*****
      SUBROUTINE VBDECP(A,LL,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(1,1),LL(1)
C**
C**   CHOLESKY DECOMPOSITION A=LDU
C**
      JM=1
      DO 40 J=1,N
      JA=LL(J)
      JM=J-JA+JM
      J1=J-1
      JM1=JM+1
      IF(JM1.GT.J1) GO TO 25
      DO 20 I=JM1,J1
      I1=I-1
      IA=LL(I)
      IM=MAX0(I-IA+LL(I1),JM)
      AIJ=0.0
      IF(IM.GT.I1) GO TO 20
      DO 10 K=IM,I1
   10 AIJ=AIJ+A(K,IA)*A(K,JA)
   20 A(I,JA)=A(I,JA)-AIJ
   25 IF(JM.GT.J1) GO TO 35
      AJJ=0.0
      DO 30 I=JM,J1
      IA=LL(I)
      AIJ=A(I,JA)/A(I,IA)
      AJJ=AJJ+AIJ*A(I,JA)
   30 A(I,JA)=AIJ
      A(J,JA)=A(J,JA)-AJJ
   35 IF(A(J,JA).EQ.0.0) GO TO 50
   40 JM=JA
      RETURN
   50 PRINT 85,J,J
   85 FORMAT(3H A( I3,1H,I3,5H)=0.0/)
      RETURN
      END

      SUBROUTINE VBSOLX(A,B,LL,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(1,1),B(1),LL(1)
C**
C**   AX=B  :  LDUX=B  ;  WHERE L=U'
C**   FORWARD SUBSTITUTION LZ=B
C**
      IM=1
      DO 70 I=1,N
      I1=I-1
      IA=LL(I)
      IM=I-IA+IM
      IF(IM.GT.I1) GO TO 70
      BI=B(I)
      DO 60 K=IM,I1
   60 BI=BI-A(K,IA)*B(K)
      B(I)=BI
   70 IM=IA
C**
C**   DIVIDED BY DIAGONAL DY=Z
C**
      DO 75 I=1,N
      IA=LL(I)
   75 B(I)=B(I)/A(I,IA)
C**
C**   BACKWARD SUBSTITUTION UX=Y
C**
      N2=N+2
      DO 90 IN=2,N
      I=N2-IN
      I1=I-1
      IA=LL(I)
      IM=I-IA+LL(I1)
      IF(IM.GT.I1) GO TO 90
      BI=B(I)
      DO 80 K=IM,I1
   80 B(K)=B(K)-A(K,IA)*BI
   90 CONTINUE
      RETURN
      END
C*****
      SUBROUTINE BVEXAV(A,B,V,N,X,IH,KK,EPSI,NV)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(1),B(1),V(1,1),X(1),IH(1)
C**    B(N,N) * V(N,J) = A(N,N) * V(N,J) * X(J)
C**   VT(J,N) * A(N,N) * V(N,J) = 1
C**   VT(J,N) * B(N,N) * V(N,J) = X(J)
      CALL JACOBI(B,V,N,X,IH,1,KK,EPSI,NV)
      DO 40 J=1,N
      JV=(J-1)*NV+1
C     JJ=(J-1)*NA+J
      JJ=(J+1)*J/2
      IF(B(JJ).LE.0.0) STOP 1
      BB=1.0/SQRT(B(JJ))
      DO 40 I=1,N
   40 V(I,JV)=V(I,JV)*BB
      CALL MAV(A,B,A,V,X,N,N,NV)
      CALL MAV(A,B,B,V,X,N,0,NV)
      CALL JACOBI(A,V,N,X,IH,-1,KK,EPSI,NV)
      DO 50 J=1,N
      JV=(J-1)*NV+1
C     JJ=(J-1)*NA+J
      JJ=(J+1)*J/2
      IF(A(JJ).EQ.0.0) STOP 2
      X(J)=1.0/A(JJ)
      BB=SQRT(ABS(X(J)))
      DO 50 I=1,N
   50 V(I,JV)=V(I,JV)*BB
      RETURN
      END
C*****
      SUBROUTINE MAV(A,B,C,V,X,N,M,NV)
C**      CALL    MAV(A,B,A,V,X,N,N,NV)
C**      CALL    MAV(A,B,B,V,X,N,0,NV)
C**   (CT+A) *    V   --> (AT+B) ; B UNUSED IF M=0
C**   A A A     V V V     A B B
C**   C A A  *  V V V  =  A A B
C**   C C A     V V V     A A A
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(1,1),B(1,1),C(1,1),V(1,1),X(1)
      IA=1
      DO 50 I=1,N
      I1=I-1
      DO 10 J=1,I1
   10 X(J)=C(J,IA)
      JA=IA
      DO 20 J=I,N
      X(J)=A(I,JA)
C  20 JA=JA+NA
   20 JA=JA+J
      JV=1
      DO 30 J=1,I
      AV=0.0
      DO 25 K=1,N
   25 AV=AV+X(K)*V(K,JV)
      A(J,IA)=AV
   30 JV=JV+NV
      IF(I.GE.M) GO TO 50
      I1=I+1
C     JB=IA+NA
      JB=IA+I
      DO 40 J=I1,N
      AV=0.0
      DO 35 K=1,N
   35 AV=AV+X(K)*V(K,JV)
      B(I,JB)=AV
C     JB=JB+NA
      JB=JB+J
   40 JV=JV+NV
C  50 IA=IA+NA
   50 IA=IA+I
      RETURN
      END
C*****
      SUBROUTINE JACOBI(A,V,N,X,IH,JEVC,M,EPSI,NV)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(1,1),V(1,1),X(1),IH(1)
      IF(JEVC) 15,15,10
   10 DO 14 J=1,N
      JV=(J-1)*NV+1
      DO 12 I=1,N
   12 V(I,JV)=0.0
   14 V(J,JV)=1.0
C**
   15 M=M
      IF(N-1) 400,400,20
   20 AMIN=ABS(A(1,1))
      DO 40 JX=2,N
      J1=JX-1
C     JA=J1*NA+1
      JA=JX*J1/2+1
      II=JX
      IF(AMIN-ABS(A(II,JA))) 30,30,25
   25 AMIN=ABS(A(II,JA))
   30 X(JX)=-1.0
      DO 40 II=1,J1
      IF(X(JX)-ABS(A(II,JA))) 35,40,40
   35 X(JX)=ABS(A(II,JA))
      IH(JX)=II
   40 CONTINUE
C**
   50 AIJ=-1.0
      DO 60 JX=2,N
      IF(AIJ-X(JX)) 55,60,60
   55 AIJ=X(JX)
      IP=IH(JX)
      JP=JX
   60 CONTINUE
C**
      IF(AIJ-AMIN*EPSI) 400,400,70
   70 M=M+1
      IA=IP
      IB=JP
C     JA=(IP-1)*NA+1
C     JB=(JP-1)*NA+1
      JA=IP*(IP-1)/2+1
      JB=JP*(JP-1)/2+1
      AIJ=2.*A(IA,JB)
      AII=A(IA,JA)
      AJJ=A(IB,JB)
      TANG=AIJ/(ABS(AJJ-AII)+SQRT((AJJ-AII)**2+AIJ**2))
      IF(AII-AJJ) 75,80,80
   75 TANG=1.0/TANG
   80 COSS=1.0/(1.0+TANG*TANG)
      COSN=SQRT(COSS)
      SINE=TANG*COSN
      A(IA,JB)=0.0
      A(IA,JA)=((AJJ*TANG+AIJ)*TANG+AII)*COSS
      A(IB,JB)=((AII*TANG-AIJ)*TANG+AJJ)*COSS
C**
      IF(AMIN-ABS(A(IB,JB))) 90,90,85
   85 AMIN=ABS(A(IB,JB))
C**
   90 JX=IP
      JY=JP
      X(JX)=0.0
      X(JY)=0.0
      DO 250 JJ=1,N
      IF(JJ-IP) 100,110,120
  100 IA=JJ
  105 IB=JJ
      GO TO 150
  110 IA=IP
      GO TO 250
  120 JX=JJ
C     JA=JA+NA
      JA=JA+JJ-1
      IF(JJ-JP) 105,130,140
  130 IB=JP
      GO TO 250
  140 JY=JX
      JB=JA
  150 AIJ=A(IA,JA)
      A(IA,JA)= AIJ*COSN+A(IB,JB)*SINE
      A(IB,JB)=-AIJ*SINE+A(IB,JB)*COSN
      IF((IH(JX)-IP)*(IH(JY)-JP)) 200,160,200
  160 X(JX)=-1.0
      J1=JJ-1
      DO 180 II=1,J1
      IF(X(JX)-ABS(A(II,JA))) 170,180,180
  170 X(JX)=ABS(A(II,JA))
      IH(JX)=II
  180 CONTINUE
      IF(JJ-JP) 220,250,250
  200 IF(X(JX)-ABS(A(IA,JA))) 210,220,220
  210 X(JX)=ABS(A(IA,JA))
      IH(JX)=IA
  220 IF(X(JY)-ABS(A(IB,JB))) 230,250,250
  230 X(JY)=ABS(A(IB,JB))
      IH(JY)=IB
  250 CONTINUE
C**
      IF(JEVC) 300,50,300
  300 JU=(IP-1)*NV+1
      JV=(JP-1)*NV+1
      DO 310 I=1,N
      AIJ=V(I,JU)
      V(I,JU)= AIJ*COSN+V(I,JV)*SINE
  310 V(I,JV)=-AIJ*SINE+V(I,JV)*COSN
      GO TO 50
  400 RETURN
      END