program mtmsvdspec
c
c     Copyright (C) 1995
c
c     Michael Mann
c
c     MTM-SVD -- evolutive/fixed window  LFV spectrum estimation
c
c     v1.0
c
c     uses  original  multitaper spectral code of Jeffrey Park (c) 
c
c   mmax     - maximum number of spatial dofs desired (ie, max # of time series)
c   ntapmax  - maximum number of spectral dofs desired
c   ifmax    - maximum padded DFT length (maximum=8192)
c   nscanmax - maximum length of time series
c
c
      parameter (mmax=25,
     $           nscanmax=1200,ntapmax=3,ifmax=4096)
c
      character*19 filename
      character*1 cfnum
      complex *8 AAA(mmax,ntapmax,ifmax)
      complex *16 AA(mmax,ntapmax),UU(mmax,mmax),VV(ntapmax,ntapmax)
      real*8 ell
      real *8 S(ntapmax),sum
      real *8 partsum(ntapmax)
      real ar,ai
      real weight(mmax),a(nscanmax,mmax),anew(nscanmax,mmax)
      real datave(mmax),sd(mmax)
      real var
      common/npiprol/anpi
      common/nnaamm/iwflag
      common/work/b1(8194)
      common/taperz/ta(4100,16,2),tai(4100,16,2),dcf(4100,16,2),
     $          amu(4100,4)
      common/stap2/tas(8192,16)
      common/envel1/ell(20)
c
c
      iwflag=0
      zero = 0.0
      one = 1.0
      two = 2.0
      pi=3.14159265358979

c
c     dt    = time interval spacing in units of years
c     nscan = length of time series in units of "dt"
c     iabv  = number of time series in analysis (ie, "spatial" array)
c
c     presently set for 25 records of 1 month resolution and 100 year
c     length (** make sure array dimensions are large enough:
c          "ifmax" must be greater than nscan **)
c     so nscan = 100*12 = 1200, iabv = 25
c
c     dt = 1/12 so that time units/frequencies will be in terms
c     of annual time units
c
c     **** THESE SETTINGS ARE DATASET DEPENDENT!!! ****
c
      dt = 0.083333
      nscan = 1200
      iabv = 25
c
c     set frequency range for DFT, frequency spacing, zero
c     padding, etc.
c     ** NOTE: maximum padded DFT length is npad=8192 **
c
c     f1 = minimum frequency for analysis = 0
c     f2 = maximum frequency  (in cyc/year -- this can be varied up to
c          the nyquist sampling frequency = 0.5/dt)
c
c     ( presently set for the range 0 to 0.5 cycle year, or, periods greater
c      than tau=2 years. This is a smaller fraction of the full Nyquist
c      interval f=0 to f=6.0 cycle/year for monthly sampling)
c
      f1 = 0.0
      f2 = 0.5

c     default settings
c 
      npi = 2
      anpi=float(npi)
      nwin=3
      itype = 0
      nbeg=1
      nend = nscan
      neval = nscan
      nmove = nscan
      ninterv = nscan

444   write (6,*) 'present settings'
      write (6,*) '1) use',nwin, ' x ',npi,' pi tapers'
      write (6,*) '2) data interval (',nbeg,' to ',nend,' )'
      if (itype.eq.0) then 
        write (6,*) '3) fixed window analysis'
      else
         write (6,*) '3) evolutive analysis'
         write (6,*) '    --',nmove,' unit moving window'
         write (6,*) '    --',ninterv,' unit time step'
      endif
c      write (6,*) '-----------------'
      write (6,*) 'item to change (0=continue)'
      read (5,*) icont
      if (icont.eq.1) then
         goto 1111
      else
        if (icont.eq.2) then
           goto 2222
        else
           if (icont.eq.3) then
              goto 3333
           else
              goto 9999
           endif
        endif
      endif


1111  write (6,*) 'time-frequency bandwidth product "p"'
      read (5,*) npinew
      if (ntapmax.lt.2*npinew-1) then
         write (6,*) 'sorry--exceeds maximum bandwidth'
      else
         write (6,*) 'number of tapers'
         read (5,*) nwinnew
         if (nwinnew.gt.2*npinew-1) then
           write (6,*) 'sorry--cannot exceed 2*p-1'
         else
           npi=npinew
           anpi = float(npi)
           nwin=nwinnew
         endif
      endif
      goto 444

2222  write (6,*) 'beginning time index'
      read (5,*) nbeg
      write (6,*) 'ending time index'
      read (5,*) nend
      nmove = nend-nbeg+1
      neval = nmove
      ninterv = nmove
      goto 444

3333  write (6,*) 'choose'
      write (6,*) '(0) fixed window'
      write (6,*) '(1) evolutive analysis'
      read (5,*) itype
      if (itype.eq.1) then
         write (6,*) 'window width (in samples--e.g., # months)'
         read (5,*) nmove
         write (6,*) 'interval between evaluations (in samples " )'
         read (5,*) ninterv
      endif
      goto 444
c
c     npad = zero padding (power of two--must be > nscan)
c
c
9999  npad=4096
c
c     calculate frequency spacing, etc.
c
      ddf=1./(npad*dt)
      eps=.1*ddf
      nf1=(f1+eps)/ddf
      nf2=(f2+eps)/ddf
      nf=nf2-nf1+1
      nfmax = nf-1
c
c     set spatial, spectral dimensions for SVD matrix
c
      M = iabv
      N = nwin
      NV = N
      NU = N
      IP = 0
c
c     read in data.
c
c     * THIS SECTION WILL HAVE TO MODIFIED DEPENDING ON THE DATA
c       ANALYZED *
c
      open (unit=9,file='synth1.dat',
     $         status='old')
c   
      open (unit=10,file='synth2.dat',
     $         status='old')
c
      open (unit=11,file='synth3.dat',
     $         status='old')
c
      open (unit=12,file='synth4.dat',
     $         status='old')
c
      open (unit=13,file='synth5.dat',
     $         status='old')
c
c    read in the iabv=25 gridpoints from the
c    five files of 5 gridpoint series/file
c
      do j=1,nscan
          read (9,*) adum,(a(j,jnode),jnode=1,5)
          read (10,*) adum,(a(j,jnode),jnode=6,10)
          read (11,*) adum,(a(j,jnode),jnode=11,15)
          read (12,*) adum,(a(j,jnode),jnode=16,20)
          read (13,*) adum,(a(j,jnode),jnode=21,25)
      end do
      close (unit=9)
      close (unit=10)
      close (unit=11)
      close (unit=12)
      close (unit=13)
c
c     set uniform weights on each of the gridpoints. This is where
c     we might normally want to set gridpoint dependent weights
c     e.g., cos(latitude) weighting to provide areal representivity
c     to the decomposition, etc. For simplicity, we set the weights
c     as uniform here.
c
      do j=1,iabv
         weight(j)=one
      end do
c
777   format (11f6.2)
c
c
c
      write (6,*) 'read in: '
      write (6,*) iabv, ' time series'
      write (6,*) nscan, ' sequential time snapshots of data...'
c
c     for each time series, subtract mean, normalize by standard 
c     deviation for each month,
c     also set weights for the different dataseries. 
c     ( for the the present case, all dataseries are equally weighted for
c       analysis)
c
c
      open (unit=15,file='svd-means.out',
     $     status='unknown')
      open (unit=16,file='svd-sigmas.out',
     $     status='unknown')
      do j=1,iabv
         sum=0.0
         var = 0.0
         do i=nbeg,nend
            sum=sum+a(i,j)
         end do
         datave(j)=sum/float(neval)
         inew = 0
         do i=nbeg,nend
           inew = inew+1
           anew(inew,j)=a(i,j)-datave(j)
           var =var + anew(inew,j)**2
         end do
         sd(j)=sqrt(var/float(neval))
         do i=1,neval
           anew(i,j)=weight(j)*anew(i,j)/sd(j)
         end do
         write (15,*) j,datave(j)
         write (16,*) j,sd(j)
      end do 
      close (unit=15)
      close (unit=16)
c            
c     open output files
c
      do i=1,nwin
        cfnum = char(48+i)
        filename = 'evals-'//cfnum//'.out'
        open (unit=13+i-1,file=filename,
     $     status='unknown')
      end do
c
c     construct tapers appropriate for window width
c
      npts = nmove
      call taper(npts,nwin,ell)
c
      write (6,*) 'nscan ',nscan
      write (6,*) 'neval ',neval
      write (6,*) 'nbeg ',nbeg
      write (6,*) 'nend ',nend
      write (6,*) 'dt    ',dt
      write (6,*) 'nmove ',nmove
      write (6,*) 'ninterv ',ninterv
      niter = (neval-nmove)/ninterv + 1
      write (6,*) 'beginning 1st of ',niter, ' iteration/s...'
c
c     begin moving time window loop
c
      do it = 1,niter
        n1 = 1+(it-1)*ninterv
        n2 = n1+npts-1
        ncent = (n1+n2-1)/2
c
c       calculate DFTs for each gridpoint
c
        do ii=1,iabv
c        demean gridpoint series locally
         sum = 0.
         do i=n1,n2
            sum = sum+anew(i,ii)
         end do
         sum = sum/float(n2-n1+1)
         do iwin=1,nwin
          do i=1,npad
             b1(i)=0.
          end do
          do i=n1,n2
            b1(i-n1+1)=(anew(i,ii)-sum)*tas(i+1-n1,iwin)
          end do
          call realft(b1,npad/2,1)
          b1(npad-1)=b1(2)
          b1(npad)=0.
          b1(2)=0.
          do iif=1,nfmax
           iwhere = 2*(nf1+iif)-1
           ar=b1(iwhere)
           ai=b1(iwhere+1)
           AAA(ii,iwin,iif)=cmplx(ar,ai)
          end do
         end do
c
c       terminate DFT calculation
c
        end do
c
c       perform SVD for given time and frequency
c
        do iif =1,nfmax
         do ii=1,iabv
           do iwin=1,nwin
              AA(ii,iwin)=1.d0*AAA(ii,iwin,iif)
           end do
         end do
c     
c        perform SVD 
c
         call CSVD(AA,mmax,nwin,M,N,IP,NU,NV,S,UU,VV)
c
c
c        determine local fractional variances 
c
         sum = 0.d0
         do i=1,nwin
          sum = sum + S(i)**2
          partsum(i)=0.0
          do j=i,nwin
            partsum(i)=partsum(i)+S(j)**2
          end do
         end do
         freq = (iif-1)*ddf+f1
         singval = real(S(1))
         fracvar = real(S(1)**2/sum)
         if (itype.eq.1) then
          write (13,6666) freq,float(ncent+nbeg-1)*dt,
     $       fracvar,singval,iif
         else
          write (13,7777) freq,fracvar,singval,iif
         endif
         do i=2,nwin
          singval = real(S(i))
          fracvar = real(S(i)**2/sum)
          partvar = real(S(i)**2/partsum(i))
          if (itype.eq.1) then
               write (13+i-1,888) freq,float(ncent+nbeg-1)*dt,
     $           partvar,fracvar,singval,iif
          else
            write (13+i-1,999) freq, partvar,fracvar,singval,iif
          endif
         end do
c
c       terminate frequency loop
c
        end do
c  
c       provide progress update
c
        open (unit=12,file='svdevals-synth.status',
     $         status='unknown')
        if ((niter.ge.10).and.(mod(it-1,niter/10).eq.0)) then
         write (12,*) it,' iterations of ',niter,' done...'
         write (6,*) it,' iterations of ',niter,' done...'
        endif
        close (unit=12)
c
c       terminate time loop
c
      end do
      write (6,*) 'done.'
c
c     done -- close files
c
      do i=1,nwin
         close (unit=13+i-1)
      end do 
      close (unit=33)
888   format (f9.4,f6.2,2f9.4,f12.4,i6)
999   format (f9.4,2f9.4,f12.4,i6)
6666  format (f9.4,f6.2,f9.4,f12.4,i6)
7777  format (f9.4,f9.4,f12.4,i6)
      stop
      end
c
c
c
      subroutine taper(n,nwin,el)
c
c  to generate slepian tapers
c  ta is a real*4 array
c
c         j. park
c
      real*8 el,a,z,pi,ww,cs,ai,an,eps,rlu,rlb
      real*8 dfac,drat,gamma,bh,ell
      common/nnaamm/iwflag
      common/npiprol/anpi
      common/tapsum/tapsum(20),ell(20)
      common/work/ip(8194)        
      common/taperzz/z(65536)  ! we use this common block for scratch space
      common/stap2/ta(8192,16)
      dimension a(8192,8),el(10)
      data iwflag/0/,pi/3.14159265358979d0/
      equivalence (a(1,1),ta(1,1))
      an=dfloat(n)
      ww=dble(anpi)/an
      cs=dcos(2.d0*pi*ww)
c initialize matrix for eispack subroutine
c      type *,'ww,cs,an',ww,cs,an
      do i=0,n-1
        ai=dfloat(i)
        a(i+1,1)=-cs*((an-1.d0)/2.d0-ai)**2
        a(i+1,2)=-ai*(an-ai)/2.d0
        a(i+1,3)=a(i+1,2)**2        ! required by eispack routine
      end do
      eps=1.e-13
      m11=1
      call tridib(n,eps,a(1,1),a(1,2),a(1,3),rlb,rlu,m11,nwin,el,ip,
     x       ierr,a(1,4),a(1,5))
c      type *,ierr,rlb,rlu
c      call tnou('eigenvalues for the eigentapers')
c      type *,(el(i),i=1,nwin)
      call tinvit(n,n,a(1,1),a(1,2),a(1,3),nwin,el,ip,z,ierr,
     x            a(1,4),a(1,5),a(1,6),a(1,7),a(1,8))      
c  we calculate the eigenvalues of the dirichlet-kernel problem
c  i.e. the bandwidth retention factors
c  from slepian 1978 asymptotic formula, gotten from thomson 1982 eq 2.5
c  supplemented by the asymptotic formula for k near 2n from slepian 1978 eq 61
      dfac=an*pi*ww
      drat=8.d0*dfac
      dfac=4.d0*dsqrt(pi*dfac)*dexp(-2.d0*dfac)
      do k=1,nwin
        el(k)=1.d0-dfac
        dfac=dfac*drat/k  ! is this correct formula? yes,but fails as k -> 2n
      end do
c      call tnou('eigenvalues for the eigentapers (small k)')
c      type *,(el(i),i=1,nwin)
      gamma=dlog(8.d0*an*dsin(2.d0*pi*ww))+0.5772156649d0
      do k=1,nwin
        bh=-2.d0*pi*(an*ww-dfloat(k-1)/2.d0-.25d0)/gamma
        ell(k)=1.d0/(1.d0+dexp(pi*bh))
      end do
c      call tnou('eigenvalues for the eigentapers (k near 2n)')
c      type *,(ell(i),i=1,nwin)
      do i=1,nwin
        el(i)=dmax1(ell(i),el(i))
      end do     
c      call tnou('composite asymptotics for taper eigenvalues')
c      type *,(el(i),i=1,nwin)  
c   normalize the eigentapers to preserve power for a white process
c   i.e. they have rms value unity
      do k=1,nwin
        kk=(k-1)*n
        tapsum(k)=0.
        tapsq=0.
        do i=1,n
          aa=z(kk+i)
          ta(i,k)=aa
          tapsum(k)=tapsum(k)+aa
          tapsq=tapsq+aa*aa
        end do
        aa=sqrt(tapsq/n)
        tapsum(k)=tapsum(k)/aa
        do i=1,n
          ta(i,k)=ta(i,k)/aa
        end do
      end do
c      type *,'tapsum',(tapsum(i),i=1,nwin)
c  refft preserves amplitudes with zeropadding
c  zum beispiel: a(i)=1.,i=1,100 will transform at zero freq to b(f=0.)=100
c  no matter how much zero padding is done
c  therefore we need not doctor the taper normalization,
c  but wait until the fft to force the desired units
      iwflag=1
      return
      end
c
c
c
      SUBROUTINE CSVD (A,MMAX,NMAX,M,N,IP,NU,NV,S,U,V)
C
C   Singular value decomposition of an  M by N  complex matrix  A,
C   where M .GT. N .  The singular values are stored in the vector
C   S. The first NU columns of the M by M unitary matrix U and the
C   first NV columns of the N by N unitary matrix V  that minimize
C   Det(A-USV*) are also computed.
C
C
C         P.A. Businger and G.H. Golub, "Singular Value Decomposition
C         of a Complex Matrix," Communications of the ACM, vol. 12,
C         pp. 564-565, October 1969.
C
C   This algorithm is reprinted by permission, Association for
C   Computing Machinery; copyright 1969.
C
      parameter (nwork=10000)
      COMPLEX *16 A(MMAX,NMAX),U(MMAX,MMAX),V(NMAX,NMAX),Q,R
      REAL *8 S(NMAX),B(nwork),C(nwork),T(nwork),zero,one,ETA,TOL,
     $        EPS
      ETA = 1.2d-7
      TOL = 2.4d-32
      zero = 0.d0
      one = 1.d0
      NP=N+IP
      N1=N+1
C   Householder reduction
      C(1)=zero
      K=1
10    K1=K+1
C   Elimination of A(I,K) , I=K+1,...,M
      Z=zero
      DO 20 I=K,M
20      Z=Z+REAL(A(I,K))**2+dIMAG(A(I,K))**2
      B(K)=zero
      IF (Z .LE. TOL)  GO TO 70
      Z=SQRT(Z)
      B(K)=Z
      W=CdABS(A(K,K))
      Q=dcmplx(one,zero)
      IF (W .NE. zero)  Q=A(K,K)/W
      A(K,K)=Q*(Z+W)
      IF (K .EQ. NP)  GO TO 70
      DO 50 J=K1,NP
        Q=dcmplx(zero,zero)
        DO 30 I=K,M
30        Q=Q+dCONJG(A(I,K))*A(I,J)
        Q=Q/(Z*(Z+W))
        DO 40 I=K,M
40        A(I,J)=A(I,J)-Q*A(I,K)
50      CONTINUE
C   Phase transformation
      Q=-dCONJG(A(K,K))/CdABS(A(K,K))
      DO 60 J=K1,NP
60      A(K,J)=Q*A(K,J)
C   Elimination of A(K,J) , J=K+2,...,N
70    IF (K .EQ. N)  GO TO 140
      Z=zero
      DO 80 J=K1,N
80      Z=Z+REAL(A(K,J))**2+dIMAG(A(K,J))**2
      C(K1)=zero
      IF (Z .LE. TOL)  GO TO 130
      Z=SQRT(Z)
      C(K1)=Z
      W=CdABS(A(K,K1))
      Q=dcmplx(one,zero)
      IF (W .NE. zero)  Q=A(K,K1)/W
      A(K,K1)=Q*(Z+W)
      DO 110 I=K1,M
        Q=dcmplx(zero,zero)
        DO 90 J=K1,N
90        Q=Q+dCONJG(A(K,J))*A(I,J)
        Q=Q/(Z*(Z+W))
        DO 100 J=K1,N
100       A(I,J)=A(I,J)-Q*A(K,J)
110     CONTINUE
C   Phase transformation
      Q=-dCONJG(A(K,K1))/CdABS(A(K,K1))
      DO 120 I=K1,M
120     A(I,K1)=A(I,K1)*Q
130   K=K1
      GO TO 10
C   Tolerance for negligible elements
140   EPS=zero
      DO 150 K=1,N
        S(K)=B(K)
        T(K)=C(K)
150   EPS=DMAX1(EPS,S(K)+T(K))
      EPS=EPS*ETA
C   Initialization of U and V
      IF (NU .EQ. 0)  GO TO 180
      DO 170 J=1,NU
        DO 160 I=1,M
160       U(I,J)=dcmplx(zero,zero)
170     U(J,J)=dcmplx(one,zero)
180   IF (NV .EQ. 0)  GO TO 210
      DO 200 J=1,NV
        DO 190 I=1,N
190       V(I,J)=dcmplx(zero,zero)
200     V(J,J)=dcmplx(one,zero)
C   QR diagonalization
210   DO 380 KK=1,N
        K=N1-KK
C   Test for split
220     DO 230 LL=1,K
          L=K+1-LL
          IF (ABS(T(L)) .LE. EPS)  GO TO 290
          IF (ABS(S(L-1)) .LE. EPS)  GO TO 240
230       CONTINUE
C   Cancellation of B(L)
240     CS=zero
        SN=one
        L1=L-1
        DO 280 I=L,K
          F=SN*T(I)
          T(I)=CS*T(I)
          IF (ABS(F) .LE. EPS)  GO TO 290
          H=S(I)
          W=SQRT(F*F+H*H)
          S(I)=W
          CS=H/W
          SN=-F/W
          IF (NU .EQ. 0)  GO TO 260
          DO 250 J=1,N
            X=REAL(U(J,L1))
            Y=REAL(U(J,I))
            U(J,L1)=dCMPLX(X*CS+Y*SN,0.)
250         U(J,I)=dCMPLX(Y*CS-X*SN,0.)
260       IF (NP .EQ. N)  GO TO 280
          DO 270 J=N1,NP
            Q=A(L1,J)
            R=A(I,J)
            A(L1,J)=Q*CS+R*SN
270         A(I,J)=R*CS-Q*SN
280       CONTINUE
C   Test for convergence
290     W=S(K)
        IF (L .EQ. K)  GO TO 360
C   Origin shift
        X=S(L)
        Y=S(K-1)
        G=T(K-1)
        H=T(K)
        F=((Y-W)*(Y+W)+(G-H)*(G+H))/(2.d0*H*Y)
        G=SQRT(F*F+one)
        IF (F .LT. zero)  G=-G
        F=((X-W)*(X+W)+(Y/(F+G)-H)*H)/X
C   QR step
        CS=one
        SN=one
        L1=L+1
        DO 350 I=L1,K
          G=T(I)
          Y=S(I)
          H=SN*G
          G=CS*G
          W=SQRT(H*H+F*F)
          T(I-1)=W
          CS=F/W
          SN=H/W
          F=X*CS+G*SN
          G=G*CS-X*SN
          H=Y*SN
          Y=Y*CS
          IF (NV .EQ. 0)  GO TO 310
          DO 300 J=1,N
            X=REAL(V(J,I-1))
            W=REAL(V(J,I))
            V(J,I-1)=dCMPLX(X*CS+W*SN,0.)
300         V(J,I)=dCMPLX(W*CS-X*SN,0.)
310       W=SQRT(H*H+F*F)
          S(I-1)=W
          CS=F/W
          SN=H/W
          F=CS*G+SN*Y
          X=CS*Y-SN*G
          IF (NU .EQ. 0)  GO TO 330
          DO 320 J=1,N
            Y=REAL(U(J,I-1))
            W=REAL(U(J,I))
            U(J,I-1)=dCMPLX(Y*CS+W*SN,0.)
320         U(J,I)=dCMPLX(W*CS-Y*SN,0.)
330       IF (N .EQ. NP)  GO TO 350
          DO 340 J=N1,NP
            Q=A(I-1,J)
            R=A(I,J)
            A(I-1,J)=Q*CS+R*SN
340         A(I,J)=R*CS-Q*SN
350       CONTINUE
        T(L)=zero
        T(K)=F
        S(K)=X
        GO TO 220
C   Convergence
360     IF (W .GE. zero)  GO TO 380
        S(K)=-W
        IF (NV .EQ. 0)  GO TO 380
        DO 370 J=1,N
370       V(J,K)=-V(J,K)
380     CONTINUE
C   Sort singular values
      DO 450 K=1,N
        G=-one
        J=K
        DO 390 I=K,N
          IF (S(I) .LE. G)  GO TO 390
          G=S(I)
          J=I
390       CONTINUE
        IF (J .EQ. K)  GO TO 450
        S(J)=S(K)
        S(K)=G
        IF (NV .EQ. 0)  GO TO 410
        DO 400 I=1,N
          Q=V(I,J)
          V(I,J)=V(I,K)
400       V(I,K)=Q
410     IF (NU .EQ. 0)  GO TO 430
        DO 420 I=1,N
          Q=U(I,J)
          U(I,J)=U(I,K)
420       U(I,K)=Q
430     IF (N .EQ. NP)  GO TO 450
        DO 440 I=N1,NP
          Q=A(J,I)
          A(J,I)=A(K,I)
440       A(K,I)=Q
450     CONTINUE
C   Back transformation
      IF (NU .EQ. 0)  GO TO 510
      DO 500 KK=1,N
        K=N1-KK
        IF (B(K) .EQ. zero)  GO TO 500
        Q=-A(K,K)/CdABS(A(K,K))
        DO 460 J=1,NU
460       U(K,J)=Q*U(K,J)
        DO 490 J=1,NU
          Q=dcmplx(zero,zero)
          DO 470 I=K,M
470         Q=Q+dCONJG(A(I,K))*U(I,J)
          Q=Q/(CdABS(A(K,K))*B(K))
          DO 480 I=K,M
480         U(I,J)=U(I,J)-Q*A(I,K)
490       CONTINUE
500     CONTINUE
510   IF (NV .EQ. 0)  GO TO 570
      IF (N .LT. 2)  GO TO 570
      DO 560 KK=2,N
        K=N1-KK
        K1=K+1
        IF (C(K1) .EQ. zero)  GO TO 560
        Q=-dCONJG(A(K,K1))/CdABS(A(K,K1))
        DO 520 J=1,NV
520       V(K1,J)=Q*V(K1,J)
        DO 550 J=1,NV
          Q=dcmplx(zero,zero)
          DO 530 I=K1,N
530         Q=Q+A(K,I)*V(I,J)
          Q=Q/(CdABS(A(K,K1))*C(K1))
          DO 540 I=K1,N
540         V(I,J)=V(I,J)-Q*dCONJG(A(K,I))
550       CONTINUE
560     CONTINUE
570   RETURN
      END