C***********************************************************************
C*    main                                                  *
C***********************************************************************
C  Source Spectrum Fixed 
      IMPLICIT REAL*8 (A-H,O-Z)
      real mm
      DIMENSION x(2583),W(2583),rat(17),y(2583),q1(17)
      DIMENSION mm(16),F1(38318),q(38318),bc(4508)
C
c 1hz-10hz 
      parameter (nm=114971,km=38335)
      dimension u(km),v(km),am(nm),ab(km),ac(nm)
      common/main/nrow
      common/main2/am,ab,ac
        m=38335
        n=2583

        OPEN(15, file='moment1.out16')
        OPEN(16, file='site11.out16')
        OPEN(17, file='qout16')
 
        OPEN(11, file='16data1.inp',status='old')
        OPEN(12, file='dmat16.inp',status='old')
        OPEN(13, file='qdata16.inp',status='old')
        OPEN(14, file='tr10r.inp',status='old')

        READ(11,*)(F1(i),i=1,38318)
        READ(12,*)(bc(j),j=1,4508)
        READ(13,*)(q(k),k=1,38318)
        READ(14,*)(rat(l),l=1,17)

         wt1=  10000
c         wt1=  1.d0
         do 1 l=1,38318
 1           u(l)=F1(l)

            do 5 i=1,17
                u(38318+i)=log10(rat(i)) * wt1
 5            continue
           
c  position of non zero in coefficient matrix
           do 11 i=1,17
            do 12 j=1,2254
               ac((3*j-2)+6762*(i-1))=bc(2*j-1)
               ac((3*j-1)+6762*(i-1))=bc(2*j)+150*(i-1)
 12            continue
 11            continue

          do 15 k=1,17
            do 16 l=1,2254
                      ac(3*l+6762*(k-1))=2566+k
 16               continue
 15           continue
 
              do 31 k=1,17
 31                   ac(114954+k)=26+150*(k-1)
c                      write(6,*)(ac(114954+k),k=1,17)
 

c  coefficient in matrix
              do 41 k=1,17
                do 42 i=1,2254
                      am((3*i-2)+6762*(k-1))=1
                      am((3*i-1)+6762*(k-1))=1
 42                   continue
 41                   continue

                do 45 i=1,17
                      do 46 j=1,2254
                      am(3*j+6762*(i-1))=-q(j+2254*(i-1))
 46                   continue
 45               continue

           do 9514 i=1,200000
           write(78,*)am(i)
9514       continue
           close(78)




               do 51 k=1,17
 51                   am(114954+k)= 1 * wt1


c  number of non zero in coefficient matrix 
c             do 80 i=1,19
c               do 81 j=1,2254
c                ab(2254*(i-1)+j)=3
c 81        continue
c 80        continue
              do 82 i=1,38318
 82               ab(i)=3
           do 85 i=1,17
                ab(38318+i)=1
 85         continue
           
          call pstomo(m,n,x,u,v,w)
       DO 100 I=1,16
           mm(I)=10**(x(I))

 100         CONTINUE
c           write(6,*)(rat(i),i=1,19)
           do 110 l=1,2550 
              y(l)=10**(x(l+16))
 110          continue
             do i=1,17
                 q1(i)=1/x(i+2566)
                 enddo

              WRITE (6,1002)(mm(L),L=1,16)
       WRITE (15,1003)(x(k),k=1,16)      
       
       write(6,1000)(y(10+(i-1)*150),rat(i),i=1,17)
       write(6,1003)(1/q1(i),i=1,17)

c         WRITE (15,1003)(x(k),k=1,17)
         write(16,1001)(y(l),l=1,2550)
         write(17,1003)(q1(l),l=1,17)
       STOP

 1000  FORMAT (1h,2f10.5) 
 1001  format(150f10.5)
 1002  format(4e10.5)
 1003  format(1f10.4)
      END


c...........................................................................
      subroutine pstomo(m,n,x,u,v,w)
      IMPLICIT REAL*8 (A-H,O-Z)
      dimension x(n),u(m),v(m),w(n)
      ITMAX = M+N+50
      RELPR = 1.E-16
      ATOL = RELPR
      BTOL = ATOL
      CONLIM = 1./RELPR
      ZERO   = 0.0
      ONE    = 1.0
      CTOL   = ZERO
      IF (CONLIM .GT. ZERO) CTOL = ONE/CONLIM
      DAMPSQ = DAMP**2
      ANORM  = ZERO
      ACOND  = ZERO
      BBNORM = ZERO
      DDNORM = ZERO
      RES2   = ZERO
      XNORM  = ZERO
      XXNORM = ZERO
      CS2    = -ONE
      SN2    = ZERO
      Z      = ZERO
      ITN    = 0
      ISTOP  = 0
      NSTOP  = 0
c      
      do i=1,n 
        x(i)=0.
        v(i)=0. 
      enddo 
      call normlz(m,u,beta)
      b1=beta
      call atupv(m,n,u,v)
      call normlz(n,v,alfa)
      rhobar=alfa
      phibar=beta
c
      BNORM  = BETA
      RNORM  = BETA
      ARNORM = ALFA*BETA
      IF (ARNORM .LE. ZERO) GO TO 800
      TEST1  = ONE
      TEST2  = ALFA/BETA
c
      do i=1,n 
        w(i)=v(i) 
      enddo
      do 99 iter =1,itmax 
        a=-alfa
        do i=1,m 
          u(i)=a*u(i) 
        enddo
        call avpu(m,n,u,v)
        call normlz(m,u,beta)
c        BBNORM = BBNORM + ALFA**2 + BETA**2 + DAMPSQ
        BBNORM = BBNORM + ALFA**2 + BETA**2
        b=-beta
        do i=1,n 
          v(i)=b*v(i) 
        enddo
        call atupv(m,n,u,v)
        call normlz(n,v,alfa)
c        RHBAR2 = RHOBAR**2 + DAMPSQ 
        RHBAR2 = RHOBAR**2
        RHBAR1 = SQRT(RHBAR2)
        CS1    = RHOBAR/RHBAR1
c        SN1    = DAMP/RHBAR1
c        PSI    = SN1*PHIBAR
        SN1 = 0.0
        PSI = 0.0
        rho=sqrt(rhobar*rhobar+beta*beta)
        c=rhobar/rho
        s=beta/rho
        teta=s*alfa
        rhobar=-c*alfa
        phi=c*phibar
        phibar=s*phibar
        CS = C
        SN = S
        THETA = TETA
        TAU = SN*PHI
        t1=phi/rho
        t2=-teta/rho
        T3 = ONE/RHO
        do i=1,n 
          T = W(I)
          x(i)=t1*w(i)+x(i)
          w(i)=t2*w(i)+v(i)
          T = (T3*T)**2
          DDNORM=T+DDNORM
        enddo
        r=phibar/b1
        DELTA  =  SN2*RHO
        GAMBAR = -CS2*RHO
        RHS    = PHI - DELTA*Z
        ZBAR   = RHS/GAMBAR
        XNORM  = SQRT(XXNORM + ZBAR**2)
        GAMMA  = SQRT(GAMBAR**2 + THETA**2)
        CS2    = GAMBAR/GAMMA
        SN2    = THETA/GAMMA
        Z      = RHS/GAMMA
        XXNORM = XXNORM + Z**2
        ANORM  = SQRT(BBNORM)
        ACOND  = ANORM*SQRT(DDNORM)
        RES1   = PHIBAR**2
        RES2   = RES2 + PSI**2
        RNORM  = SQRT(RES1 + RES2)
        ARNORM = ALFA*ABS(TAU)
        TEST1  = RNORM/BNORM
        TEST2  = ARNORM/(ANORM*RNORM)
        TEST3  = ONE/ACOND
        T1     = TEST1/(ONE + ANORM*XNORM/BNORM)
        RTOL   = BTOL +  ATOL*ANORM*XNORM/BNORM
        T3 = ONE + TEST3
        T2 = ONE + TEST2
        T1 = ONE + T1
        IF (ITER.GE. ITMAX ) ISTOP = 7
        IF (T3  .LE. ONE   ) ISTOP = 6
        IF (T2  .LE. ONE   ) ISTOP = 5
        IF (T1  .LE. ONE   ) ISTOP = 4
        IF (TEST3 .LE. CTOL) ISTOP = 3
        IF (TEST2 .LE. ATOL) ISTOP = 2
        IF (TEST1 .LE. RTOL) ISTOP = 1
        IF (ISTOP.NE.0) GO TO 800 
   99 continue
  800 WRITE(*, 1900) ITER,ISTOP
      IF (ISTOP .EQ. 0) WRITE(*, 2000)
      IF (ISTOP .EQ. 1) WRITE(*, 2100)
      IF (ISTOP .EQ. 2) WRITE(*, 2200)
      IF (ISTOP .EQ. 3) WRITE(*, 2300)
      IF (ISTOP .EQ. 4) WRITE(*, 2400)
      IF (ISTOP .EQ. 5) WRITE(*, 2500)
      IF (ISTOP .EQ. 6) WRITE(*, 2600)
      IF (ISTOP .EQ. 7) WRITE(*, 2700)
 1900 FORMAT(/ 20H NO. OF ITERATIONS =, I6,
     1   8X, 21H STOPPING CONDITION =, I3)
 2000 FORMAT(/ 52H THE EXACT SOLUTION IS  X = 0.                      )
 2100 FORMAT(/ 52H A*X - B  IS SMALL ENOUGH, GIVEN  ATOL, BTOL        )
 2200 FORMAT(/ 52H THE LEAST-SQRS SOLN IS GOOD ENOUGH, GIVEN  ATOL    )
 2300 FORMAT(/ 52H THE ESTIMATE OF  COND(ABAR)  HAS EXCEEDED  CONLIM  )
 2400 FORMAT(/ 52H A*X - B   IS SMALL ENOUGH FOR THIS MACHINE         )
 2500 FORMAT(/ 52H THE LEAST-SQRS SOLN IS GOOD ENOUGH FOR THIS MACHINE)
 2600 FORMAT(/ 52H COND(ABAR)   SEEMS TO BE TOO LARGE FOR THIS MACHINE)
 2700 FORMAT(/ 52H THE ITERATION LIMIT HAS BEEN REACHED               )
      return
      end
c.........................................................................
      subroutine normlz(n,x,s)
      IMPLICIT REAL*8 (A-H,O-Z)
      parameter (nm=114971)
      dimension x(nm)
      s=0.
      do i=1,n 
        s=s+x(i)*x(i) 
      enddo
      s=sqrt(s)
      ss=1./s
      do i=1,n 
        x(i)=x(i)*ss 
      enddo
      return
      end
c.........................................................................
      subroutine atupv(m,n,u,v)
      IMPLICIT REAL*8 (A-H,O-Z)
      parameter (nm=114971,km=38335,nkol=2583)
      dimension u(km),v(km),am(nm),ab(km),ac(nm)
      common/main/nrow
      common/main2/am,ab,ac
      kount=0
      do 2 i=1,m 
        nc=ab(i)
        if(nc.eq.0) goto 2
        do j=1,nc 
          kount=kount+1
          jj=ac(kount)
          v(jj)=v(jj)+am(kount)*u(i)
        enddo
    2 continue
      return
      end
c........................................................................
      subroutine avpu(m,n,u,v)
      IMPLICIT REAL*8 (A-H,O-Z)
      parameter (nm=114971,km=38335)
      dimension u(km),v(km),am(nm),ab(km),ac(nm)
      common/main/nrow
      common/main2/am,ab,ac
      kount=0
      do 2 i=1,m 
        nc=ab(i)
        if(nc.eq.0)goto 2
        do j=1,nc 
          kount=kount+1
          u(i)=u(i)+am(kount)*v(ac(kount))
        enddo
    2 continue
      return
      end
c.......................................................................

