Enter your search terms Submit search form
	Web pdp-10.trailing-edge.com

C	THE PROGRAM EIG1D.F4 IS A SAMPLE PROGRAM, WHEN USED IN
C	CONJUNCTION WITH EIG1(LOADED AND SAVED AS EIG1D.SAV)
C	WILL PRODUCE A RUNNING EXAMPLE OF THE USE OF EIG1.
      SUBROUTINE EIG1(A,B,N,EPS,AM,IR,NA,NB)
	DIMENSION A(500),B(500),AM(500),IR(500)
      IJSF(I,J)=N*(I-1)+J+(I*(1-I))/2
      IJDF(I,J)=NA*(J-1)+I
    1 EN=N
    3 IF(NA-1)7,5,7
    5 MF=1
      GO TO 10
    7 MF=0
   10 IF(EPS)601,11,12
   11 SE=1.E-6
      GO TO 15
   12 SE=SQRT(EPS)
   15 DO 60 I=1,N
      DO 50 J=1,N
      IJ=NB*(J-1)+I
   50 B(IJ)=0.
      AM(I)=0.
      II=NB*(I-1)+I
   60 B(II)=1.
   80 DO 90 J=2,N
      K=J-1
      DO 90 I=1,K
      IF(MF)601,82,83
   82 L=IJDF(I,J)
      GO TO 84
   83 L=IJSF(I,J)
   84 T=ABS (A(L))
      IF(T-AM(J))90,85,85
   85 AM(J)=T
      IR(J)=I
   90 CONTINUE
  100 BA=0.
      DO 125 J=2,N
      IF(AM(J)-BA)125,120,120
  120 BA=AM(J)
      L=IR(J)
      M=J
  125 CONTINUE
      IF(MF)601,126,127
  126 LM=IJDF(L,M)
      MM=LM-L+M
      LL=IJDF(L,L)
      GO TO 138
  127 LM=IJSF(L,M)
      LL=LM-M+L
      MM=IJSF(M,M)
  138 IF(EN*ABS(A(LM))-SE)590,150,150
  150 LM1=L-1
      LP1=L+1
      MM1=M-1
      MP1=M+1
  200 R=SQRT((A(LL)-A(MM))**2+4.*A(LM)*A(LM))
      T=A(LL)-A(MM)
  400 C=SQRT(.5*(1.+ABS(T)/R))
      S=-A(LM)/(R*C)
      IF(T)401,402,402
  401 S=-S
  402 SS=S*S
      CC=C*C
      SC=S*C
      P=2.*A(LM)*SC
	A(LM)=0
	AM(L)=0
	AM(M)=0
      IF(LM1)411,450,411
  411 DO 440 I=1,LM1
      IF(MF)601,412,413
  412 IL=IJDF(I,L)
      IM=IJDF(I,M)
      GO TO 414
  413 IL=IJSF(I,L)
      IM=IL-L+M
  414 TE=A(IL)
      A(IL)=A(IL)*C-S*A(IM)
      A(IM)=A(IM)*C+TE*S
      TEM=ABS(A(IL))
      IF(TEM-AM(L))430,415,415
  415 AM(L)=TEM
      IR(L)=I
  430 TEM=ABS(A(IM))
      IF(TEM-AM(M))440,435,435
  435 AM(M)=TEM
        IR(M)=I
  440 CONTINUE
  450 IF(LP1-M)451,500,451
  451 DO 490 K=LP1,MM1
      IF(MF)601,462,463
  462 LK=IJDF(L,K)
      LT=LK-L
      KM=IJDF(K,M)
      GO TO 464
  463 LK=IJSF(L,K)
      KM=IJSF(K,M)
  464 TE=A(LK)
      A(LK)=A(LK)*C-A(KM)*S
      A(KM)=A(KM)*C+TE*S
      IF(IR(K)-L)465,470,465
  465 TEM=ABS(A(LK))
      IF(TEM-AM(K))480,467,467
  467 AM(K)=TEM
      IR(K)=L
      GO TO 480
  470 AM(K)=0.
      KM1=K-1
      DO 475 IT=1,KM1
      IF(MF)601,471,472
  471 ITK=LT+IT
      GO TO 473
  472 ITK=IJSF(IT,K)
  473 TEM=ABS(A(ITK))
      IF(TEM-AM(K))475,474,474
  474 AM(K)=TEM
      IR(K)=IT
  475 CONTINUE
  480 TEM=ABS(A(KM))
      IF(TEM-AM(M))490,485,485
  485 AM(M)=TEM
      IR(M)=K
  490 CONTINUE
  500 IF(M-N)511,580,511
  511 DO 540 J=MP1,N
      IF(MF)601,512,513
  512 LJ=IJDF(L,J)
      LT=LJ-L
      MJ=LT+M
      GO TO 514
  513 LJ=IJSF(L,J)
      MJ=IJSF(M,J)
  514 TE=A(LJ)
      A(LJ)=A(LJ)*C-A(MJ)*S
      A(MJ)=A(MJ)*C+TE*S
      IF(IR(J)-L)515,525,515
  515 IF(IR(J)-M)517,525,517
  517 TEM=ABS(A(LJ))
      IF(TEM-AM(J))520,518,518
  518 AM(J)=TEM
      IR(J)=L
  520 TEM=ABS(A(MJ))
      IF(TEM-AM(J))540,521,521
  521 AM(J)=TEM
      IR(J)=M
      GO TO 540
  525 AM(J)=0.
      JM1=J-1
      DO 530 IT=1,JM1
      IF(MF)601,526,527
  526 ITJ=LT+IT
      GO TO 528
  527 ITJ=IJSF(IT,J)
  528 TEM=ABS(A(ITJ))
      IF(TEM-AM(J))530,529,529
  529 AM(J)=TEM
      IR(J)=IT
  530 CONTINUE
  540 CONTINUE
  580 TE=A(LL)
      A(LL)=A(LL)*CC-P+A(MM)*SS
      A(MM)=A(MM)*CC+P+TE*SS
  581 DO 585 I=1,N
      IL=NB*(L-1)+I
      IM=NB*(M-1)+I
      TE=B(IL)
      B(IL)=B(IL)*C-B(IM)*S
      B(IM)=B(IM)*C+TE*S
  585 CONTINUE
      GO TO 100
590	IF(MF) 601,601,602
 602  DO 600 I=2,N
      II=IJSF(I,I)
  600 A(I)=A(II)
  601 END