c     Program for Classical Multdimensional Scaling
c             aka Principal Coordinate Analysis


c     Declare the number of objects and the number of target
c     dimensions.  Subsequent declarations depend on these
c     quantities.
c       natoms is the number of objects
c       ndim   is the number of target dimensions, e.g., 2

      integer          natoms, ndim
      parameter        (natoms=10, ndim=2)


      integer          i, j, k
      double precision tau(natoms,natoms), v(natoms), avg, 
     +                 x(ndim), stress, sum, dist


c     Declare the variables needed by ARPACK.
c       A description of these variables is given in the dssimp.f
c       driver.  Note that ncv >= 2*nev is recommended, but that
c       nev < ncv <= natoms is required.  The suggested settings 
c       should work well if ndim < natoms/2, as is virtually
c       always the case in applications of MDS.

      integer          maxn, maxnev, maxncv, ldz
      parameter        (maxn=natoms, maxnev=ndim, 
     +                 maxncv=2*maxnev, ldz=natoms)

      double precision z(ldz,maxncv), workl(maxncv*(maxncv+8)),
     +                 workd(3*maxn), d(maxncv,2), resid(maxn)
      logical          select(maxncv)
      integer          iparam(11), ipntr(11)

      character        bmat*1, which*2
      integer          ido, nar, nev, ncv, lworkl, info, ierr,
     +                 ishfts, maxitr, model
      parameter        (nar=natoms, nev=maxnev, ncv=maxncv)
      logical          rvec
      double precision tol, sigma

      double precision zero, time, time1, time2
      parameter        (zero = 0.0D+0)

      bmat = 'I'
      which = 'LA'
      lworkl = ncv*(ncv+8)
      tol = zero
      info = 0
      ido = 0
      ishfts = 1

c     Specify maxitr, the maximum number of Arnoldi iterations.
      maxitr = 300

      model = 1
      iparam(1) = ishfts
      iparam(3) = maxitr
      iparam(7) = model


c     Read the proximities from the designated file.

      print *,'Reading the proximities...'
      print *
      call timer(time1)

      open (2,file='helm.dat')
      do 81 i=1,natoms
         read (2,*) (tau(i,j), j=1,natoms)
  81  continue
      close(2)
      
      call timer(time2)
      time=time2-time1
      print *,'Number of objects = ',natoms
      print *,'Target dimension  = ',ndim
      print *
      print *,'Time =',time
      print *
      print *,'=============================='
      print *
      print *,'Computing the configuration...'
      print *
      call timer(time1)


c     Calculate the squared dissimilarities.

      do 21 i=1,natoms
         tau(i,i)=zero
  21  continue
      do 22 i=2,natoms
         do 23 j=1,i-1

c     Transform the proximities as desired.
c     For example, if the proximities are similarities,
c     with diagonal entries equal to 1 and off-diagonal
c     entries in [0,1], then it is natural to define
c     dissimilarity = 1 - similarity.
c            avg=1.0-tau(i,j)
            avg=tau(i,j)

c     Square the transformed dissimilarities.
            avg=avg**2

            tau(i,j)=avg
            tau(j,i)=avg
  23     continue
  22  continue


c     Double-center the dissimilarities.

         avg=zero
         do 31 i=1,natoms
            v(i)=0.
            do 32 j=1,natoms
               v(i)=v(i)+tau(i,j)
  32        continue
            v(i)=v(i)/natoms
            avg=avg+v(i)
  31     continue
         avg=avg/natoms
         tau(1,1) = -0.5*(tau(1,1)-v(1)-v(1)+avg)
         do 33 i=2,natoms
            do 34 j=1,i-1
               tau(i,j) = -0.5*(tau(i,j)-v(i)-v(j)+avg)
               tau(j,i) = tau(i,j)
  34        continue
            tau(i,i) = -0.5*(tau(i,i)-v(i)-v(i)+avg)
  33     continue


c        Compute the spectral decomposition using ARPACK.

  40     continue
         call dsaupd(ido,bmat,nar,which,nev,tol,resid,
     +        ncv,z,ldz,iparam,ipntr,workd,workl,
     +        lworkl,info)
         if (ido .eq. -1 .or. ido .eq. 1) then
            do 41 i=1,natoms
               avg = zero
               do 42 j=1,natoms
                  k = ipntr(1)+j-1
                  avg = avg + tau(i,j)*workd(k)
  42           continue
               k = ipntr(2)+i-1
               workd(k) = avg
  41        continue
            go to 40
         endif

         if (info .lt. 0) then
            print *,'Error with dsaupd:  info = ',info
            go to 9000
         else
            rvec = .true.
            call dseupd (rvec,'All',select,d,z,ldz,sigma,
     +           bmat,nar,which,nev,tol,resid,ncv,z,ldz,
     +           iparam,ipntr,workd,workl,lworkl,ierr)
            if (ierr .ne. 0) then
               print *, 'Error with dseupd:  ierr = ', ierr
               go to 9000
            endif

c     Print the ndim largest eigenvalues (in descending order).

            print *,'After',iparam(3),' Arnoldi iterations,'
            print *,'the',ndim,' largest eigenvalues are:'
            print *
            do 43 i=ndim,1,-1
               v(i)=d(i,1)
               print *,v(i)
  43        continue

            info = 1
            ido = 0
         endif

      call timer(time2)
      time=time2-time1
      print *
      print *,'Time =',time
      print *
      print *,'=============================='
      print *
      print *,'Writing the Cartesian coordinates...'
      print *
      call timer(time1)

c     Compute/Output the configuration.
c     The following commands write the Cartesian coordinates
c        to a file that can be read by CrystalVision.

      open (2,file='helm.cv')
      write (2,9011) ndim
9011  format('variables: ',I2)
      write (2,9012)
9012  format('labels:')
      do 360 k=1,ndim
         write (2,9013) k
 360  continue
9013  format('pc',I2.2)

      do 361 i=1,natoms
         do 362 k=1,ndim
            if (v(k) .gt. 0.) then
               x(k)=dsqrt(v(k))*z(i,k)
            else
               x(k)=zero
            endif
 362     continue

c     Print the coordinates of the current object.
c       x(ndim) is the coordinate of the most important
c       dimension,...,x(1) of the least.
c     Rewrite this statement to achieve desired format.
         write (2,*) (x(k), k=ndim,1,-1)

 361  continue

      close(2)
            
      call timer(time2)
      time=time2-time1
      print *,'Time =',time
      print *


9000  continue
      print *
      stop
      end


c====================== Subroutine timer ==============================
        
      subroutine timer(ttime)
      double precision ttime
c     *********
c
c     Subroutine timer
c     
c     This subroutine is used to determine user time. In a typical
c     application, the user time for a code segment requires calls
c     to subroutine timer to determine the initial and final time.
c
c     The subroutine statement is
c
c       subroutine timer(ttime)
c
c     where
c
c       ttime is an output variable which specifies the user time.
c
c     Argonne National Laboratory and University of Minnesota.
c     MINPACK-2 Project.
c     
c     Modified October 1990 by Brett M. Averick.
c     
c     **********
      real temp
      real tarray(2)
      real etime

c     The first element of the array tarray specifies user time   

      temp = etime(tarray)
 
      ttime = dble(tarray(1))
 
      return
 
      end  
 
c====================== The end of timer ===============================