c     Extreme Multidimensional Scaling   
c     Compute initial configuration by classical MDS;
c       then decrease the unweighted raw stress criterion
c       by Guttman majorization (GMA)

c     The user must specify the following quantities:
c       natoms = number of objects
c       ndim = target dimension
c       maxitr = number of Arnoldi iterations for CMDS
c       mgutt = number of GMA iterations

c-----DECLARATIONS------------------------------------------

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)

c     Declare the variables needed to iterate the Guttman
c     transform.
c       igutt is a counter
c       mgutt is the maximum number of iterations
c     Choose mgutt based on computational expense.

      integer          igutt, mgutt
      parameter        (mgutt=20)

      integer          i, j, k
      double precision tau(natoms,natoms), v(natoms), avg, 
     +                 x(natoms,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.

      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.
c     Choose maxitr based on computational expense.
c     There is no need to compute the eigenvalues/vectors
c       exactly (as for classical MDS), as we are merely
c       computing a reasonable initial configuration.

      maxitr = 300

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

 
c-----PREPROCESSING----------------------------------

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

c-----INITIAL CONFIGURATION---------------------

      print *
      print *,'=============================='
      print *
      print *,'Computing the initial 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     Assuming that the similarity matrix has
c     diagonal entries equal to 1 and off-diagonal
c     entries in [0,1], it is natural to define
c     dissimilarity = 1 - similarity.
c            avg=1.0-tau(i,j)
            avg=tau(i,j)

c     Transform the dissimilarities as desired.
c            avg=avg

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

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


c     Double-center the squared 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:'

            do 43 i=ndim,1,-1
               print *,d(i,1)
  43        continue
            print *

            info = 1
            ido = 0
         endif


c     Compute the classical configuration.
c       x(,ndim) is the coordinate of the most important
c       dimension,...,x(,1) of the least.

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

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

c-----Prepare for Majorization---------------------------------

      print *,'=============================='
      print *
      print *,'Preparing for majorization...'
      print *
      call timer(time1)


c     Recover the dissimilarities from the double-centered
c       squared dissmilarities.

      do 101 i=1,natoms
         v(i) = tau(i,i)
         tau(i,i) = zero
 101  continue

      do 102 i=2,natoms
         do 103 j=1,i-1
            dist = v(i)+v(j)-2.0*tau(i,j)
            if (dist .gt. 0.) then 
               dist=dsqrt(dist)
            else
               dist = zero
            endif
            tau(i,j) = dist
            tau(j,i) = dist
 103     continue
 102  continue


c     Compute/output the raw stress criterion.
c     This is optional.
               
         stress = zero    
         do 151 i=2,natoms    
            do 152 j=1,i-1
               dist = zero
               do 153 k=1,ndim
                  dist = dist+(x(i,k)-x(j,k))**2
 153           continue
               dist = dsqrt(dist)
               stress = stress+(dist-tau(i,j))**2
 152        continue
 151     continue
         print *,'Initial raw stress criterion =',stress
         print *

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

c-----GUTTMAN MAJORIZATION---------------------------------

      print *,'=============================='
      print *
      print *,'Beginning',mgutt,' GMA iterations...'
      print *
      call timer(time1)

c     Begin loop iterating the Guttman transform

      do 500 igutt=1,mgutt
      
c     Initialize new configuration Z.

         do 511 i=1,natoms
            do 512 j=1,ndim
               z(i,j) = zero
 512        continue
 511     continue

c     Compute Z=B(X)X

         do 521 i=2,natoms
            do 522 j=1,i-1
               dist = zero
               do 523 k=1,ndim
                  dist = dist+(x(i,k)-x(j,k))**2
 523           continue
               if (dist .gt. zero) then
                  avg = dsqrt(1.0/dist)
                  do 524 k=1,ndim
                     dist = tau(i,j) * avg * (x(i,k)-x(j,k))
                     z(i,k) = z(i,k)+dist
                     z(j,k) = z(j,k)-dist
 524              continue
               endif
 522        continue
 521     continue

c     Divide by n to obtain X = V^+ Z

         do 541 i=1,natoms
            do 542 j=1,ndim
               x(i,j) = z(i,j)/natoms
 542        continue
 541     continue

 500  continue

c     Compute/output the raw stress criterion.

         stress = zero
         do 551 i=2,natoms
            do 552 j=1,i-1
               dist = zero
               do 553 k=1,ndim
                  dist = dist+(x(i,k)-x(j,k))**2
 553           continue
               dist = dsqrt(dist)
               stress = stress+(dist-tau(i,j))**2
 552        continue
 551     continue
         print *,'Raw stress criterion =',stress

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

c-----POSTPROCESSING------------------------------------

      print *,'=============================='
      print *   
      print *,'Postprocessing...'
      print *
      call timer(time1)

c     Compute/Output Kruskal's stress criterion.
c     This can be used to compare solutions for different
c       values of ndim.

      stress = zero
      sum = zero
      do 601 i=2,natoms
         do 602 j=1,i-1
            dist = zero
            do 603 k=1,ndim
               dist = dist+(x(i,k)-x(j,k))**2
 603        continue
            sum = sum+dist
            dist = dsqrt(dist)
            stress = stress+(dist-tau(i,j))**2
 602     continue
 601  continue
      stress = dsqrt(stress/sum)
      print *,'Kruskal stress criterion =',stress
      print *
      print *,'Writing the Cartesian coordinates...'
      print *

c     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 660 k=1,ndim
         write (2,9013) k
 660  continue
9013  format('pc',I2.2)
               
      do 661 i=1,natoms

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

 661  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 ===============================