ˆκ”Κ‰»ƒ„ƒRƒr–@‚ΜƒvƒƒOƒ‰ƒ€iƒTƒuƒ‹[ƒ`ƒ“j

SUBROUTINE GNR(BM,BK,A,V,E,EPS,MM,MK,MA,MV,N,MAX,IPRINT)

      SUBROUTINE GNR(BM,BK,A,V,E,EPS,MM,MK,MA,MV,N,MAX,IPRINT)

C

      DIMENSION BM(MM,N),BK(MK,N),A(MA,N),V(MV,N),E(N)

      DIMENSION U(50,50),UI(50,50),UTRI(50,50),W(50,50)

      dimension temp1(50)

C

      MADJ=50

C

      CALL CHOLSK(BM,U,MA,MADJ,N)

C

      CALL UTI(U,UI,MADJ,MADJ,N)

C

      DO 20 I=1,N

      DO 10 J=1,N

      UTRI(J,I)=UI(I,J)

   10 CONTINUE

   20 CONTINUE

C

      CALL PROD(BK,UI,W,MK,MADJ,MADJ,N)

      CALL PROD(UTRI,W,A,MADJ,MADJ,MA,N)

C

      CALL JACOBI(A,E,V,EPS,MAX,N,MA,IPRINT)

C

      DO 40 I=1,N

      DO 30 J=1,N

      W(I,J)=V(I,J)

   30 CONTINUE

   40 CONTINUE

C

      CALL PROD(UI,W,V,MADJ,MADJ,MV,N)

 

      do m=1,n

         temp=e(m)

         do i=1,n

            temp1(i)=v(i,m)

            enddo

            emin=1e32

            do j=m,n

               abse=abs(e(j))

               if(abse.le.emin) then

                  jmin=j

                  emin=abse

               endif

            enddo

            e(m)=e(jmin)

            do i=1,n

               v(i,m)=v(i,jmin)

            enddo

            e(jmin)=temp

            do i=1,n

               v(i,jmin)=temp1(i)

            enddo

      enddo

 

      do m=1,n

         amax=0.0

         do I=1,n

            av=ABS(v(i,m))

            IF (av.gt.amax) amax=av

         enddo

         do i=1,n

            v(i,m)=v(i,m)/amax

         enddo

 

      enddo

 

      RETURN

      END

C

C

      SUBROUTINE UTI(U,H,MU,MH,N)

C

      DIMENSION U(MU,N),H(MH,N)

C

      DO 10 I=1,N

      IF(ABS(U(I,I)).LT.1.E-40) GO TO 1000

      H(I,I)=1.E0/U(I,I)

   10 CONTINUE

      NM1=N-1

      DO 40 IB=1,NM1

      I=N-IB

      IP1=I+1

      DO 30 JB=IP1,N

      J=N+IP1-JB

      S=0.E-60

      DO 20 K=IP1,J

      S=S-U(I,K)*H(K,J)

   20 CONTINUE

      H(I,J)=S/U(I,I)

   30 CONTINUE

   40 CONTINUE

      RETURN

C

 1000 CONTINUE

      WRITE(6,100) I,I,U(I,I)

      RETURN

C

  100 FORMAT(///,20X,'***** ERROR *****',/,10X,'U(',I3,',',I3,')=',

     &1PE15.7,' IS NEARLY EQUAL TO 0, RETURNED WITH NO CALCULATION.')

      END

C

C

      SUBROUTINE PROD(A,B,C,MA,MB,MC,N)

C

      DIMENSION A(MA,N),B(MB,N),C(MC,N)

C

      DO 30 I=1,N

      DO 20 K=1,N

      C(I,K)=0.E-60

      DO 10 J=1,N

      C(I,K)=C(I,K)+A(I,J)*B(J,K)

   10 CONTINUE

   20 CONTINUE

   30 CONTINUE

      RETURN

      END

C

C

      SUBROUTINE CHOLSK(A,U,MA,MU,N)

C

      DIMENSION A(MA,N),U(MU,N)

C

      IF(A(1,1).LE.0.) GO TO 60

      U(1,1)=SQRT(A(1,1))

      IF(N.EQ.1) RETURN

      DO 10 J=2,N

                U(1,J)=A(1,J)/U(1,1)

   10 CONTINUE

C

      DO 50 I=2,N

      IM1=I-1

      D=0.E-60

      DO 20 K=1,IM1

      D=D+U(K,I)**2

   20 CONTINUE

      U(I,I)=A(I,I)-D

      IF(U(I,I).LE.0.) GO TO 60

      U(I,I)=SQRT(U(I,I))

      IF(I.EQ.N) RETURN

      IP1=I+1

      DO 40 J=IP1,N

      D=0.E-60

      DO 30 K=1,IM1

      D=D+U(K,I)*U(K,J)

   30 CONTINUE

      U(I,J)=(A(I,J)-D)/U(I,I)

   40 CONTINUE

   50 CONTINUE

      RETURN

C

   60 CONTINUE

      WRITE(6,100)

      RETURN

C

  100 FORMAT(1H,//,10X,'***** ERROR AT CHOLESKI DECOMPOSITION',/,16X,

     &'RETURNED WITH NO CALCULATION *****')

      END

C

C

      SUBROUTINE JACOBI(A,E,V,EPS,MAX,N,M,IPRINT)

C      

      DIMENSION A(M,N),V(M,N),E(N)

C

      NM1=N-1

C

      DO 15 I=1,N

      DO 10 J=1,N

      V(I,J)=0.E-60

   10 CONTINUE

      V(I,I)=1.E0

   15 CONTINUE

C

      IP=1

      IQ=2

      DO 17 I=1,NM1

      IP1=I+1

      DO 16 J=IP1,N

      IF(ABS(A(I,J)).GT.EPS) GO TO 18

   16 CONTINUE

   17 CONTINUE

      GO TO 2000

   18 CONTINUE

      DO 1000 KOUNT=1,MAX

      IF(IPRINT.NE.0) WRITE(6,100) KOUNT,((A(I,J),V(I,J),J=1,N),I=1,N)

C

      IP=1

      IQ=2

      DO 21 I=1,NM1

      IP1=I+1

      DO 20 J=IP1,N

      IF(ABS(A(I,J)).LE.ABS(A(IP,IQ))) GO TO 20

      IP=I

      IQ=J

   20 CONTINUE

   21 CONTINUE

C

      IF(IPRINT.NE.0) WRITE(6,110) IP,IQ

C

      IF(ABS(A(IP,IQ)).LT.EPS) GO TO 2000

C      

      DV=A(IP,IP)-A(IQ,IQ)

C

      IF(ABS(DV).LE.5.E-15) GO TO 22

      D1=2.E0*A(IP,IQ)/DV

C

      TH= .5E0*ATAN(D1)

      GO TO 23

   22 CONTINUE

C

      TH= .785398164E0

   23 CONTINUE

C

      IF(IPRINT.NE.0) WRITE(6,120) TH

C

      C=COS(TH)

      S=SIN(TH)

C

      DO 24 J=1,N

      AIP=A(IP,J)

      AIQ=A(IQ,J)

      A(IP,J)=C*AIP+S*AIQ

      A(IQ,J)=-S*AIP+C*AIQ

   24 CONTINUE

C

      DO 25 I=1,N

      AIP=A(I,IP)

      AIQ=A(I,IQ)

      A(I,IP)=C*AIP+S*AIQ

      A(I,IQ)=-S*AIP+C*AIQ

      VIP=V(I,IP)

      VIQ=V(I,IQ)

      V(I,IP)=C*VIP+S*VIQ

      V(I,IQ)=-S*VIP+C*VIQ

   25 CONTINUE

 1000 CONTINUE

      WRITE(6,130) MAX

 2000 CONTINUE

      DO 30 I=1,N

      E(I)=A(I,I)

   30 CONTINUE

      RETURN

  100 FORMAT(1H,5X,I5/(1H,5X,1PE15.7,5X,E15.7))

  110 FORMAT(1H,5X,'P= ',I3,2X,'Q= ',I3)

  120 FORMAT(1H,5X,1PE15.7)

  130 FORMAT(1H,//5X,'***** ITERATIONS COUNT OVER MAXIMUM BOUND',I5,

     &' *****')

      END