program map

      parameter(mt= 3,ndy=19,ny=2004,sc=0.07)

c********************************************

c  Hinich spatial ideological map program  Version 3-19-2004

c  Method uses a set of scores which depend on Euclidean distances between

c   each respondent's ideal ideological point and the candidates positions

c   in the latent ideological space.

c  Reference: The Spatial Theory of Voting by Enelow & Hinich

c   Cambridge University Press (1984)

c************************************************************

c

c  The respondent ideal points are estimated after the candidate map is

c  estimated. The ideal point subroutine uses a three level approach:

c  1. If a respondent gives the highest possible score to one candidate

c   the ideal point of the respondent is placed at the estimated position

c   of that candidate. If the respondent gives the highest score to

c   several candidates, the respondent's ideal point is placed at the

c   centroid of those candidate's estimated positions.

c

c  2. If no candidate gets the highest score but one gets the lowest

c   possible score, the ideal point is determined from the low score and

c   the candidate who gets the highest score.

c   For example suppose that the highest score possible is 100 and the

c   lowest score possible is 0.

c   Suppose that the respondent gives candidate k a 0 and candidate j a

c   score of 90. Then his ideal point is a prorated point on the line

c   segment that connects the positions of candidates j & k.

c   Only one candidate can have the highest score for any respondent in

c   this proration class.

c

c  3. Otherwise a least squares fit is made to estimate the ideal point.

c    The least squares fit is the default.

c

c=================================================

c

c  Lines for the control file MAP.CNL

c  IMPORTANT: Comments are put after the symbol *

c  It is vital that the symbol * be used to indicate the start of the comment.

c  This delimeter must be placed in a column < = 80. Do not exceed the 80

c  column field for control statements. The parameters and control character

c  statements must be to the left of *, the delimeter.

c  Each line in map.cnl must have the symbol * in each line.

c

c=====================================

c

c  Use the map-temp.cnl as a template for the control file map.cnl

c  Use the keywords followed by an = sign and the setting. Some settings are

c   character strings while others are real numbers.

c

c=====================================

c

c  The following parameters must be set in up to 11 lines in map.cnl. They can be

c   set in any order providing that the keywords are used and the values are

c   place after the equals symbol =

c

c  IMPORTANT: Place the symbols # *  in each line after these parameters

c  are set. The program searches for the symbol # to delimite the parameter

c  reads followed by the read delimier symbol * after #.

c  The program then reads the comment line 'Candidate or party list'

c  and then reads the candidate/party names in the list.

c

c=====================================

c

c  Data file name

c

c=====================================

c

c  Euclidean distance=straight or square, political space dimension =1 or 2,

c   or equivalently one or two. The program looks for the two numbers

c   or words.

c

c=====================================

c

c  Voter id=yes or no. Yes will make the program read a respondent lable

c  from column one of the data file. This column of characters must be

c  in the data file if the switch is set =yes.

c

c=====================================

c

c  Candidate valence dimension=yes or no

c  A valence dimension score is a measure of a positive candidate/party

c  characteristic such as integrity or competence to govern.

c

c   If yes is set then a values must be placed as shown in map-temp.cnl

c   after the candidate (or party) name.

c

c=====================================

c

c  Output files=yes or no. Yes will open files Choices.graph & Voters.graph

c   which are data files to be used by any graphics package.

c

c  Choices.graph has the candidate (or party) positions estimated by map.

c  Voters.graph has the estimated voter ideal points.

c

c=====================================

c

c  Set mean=yes to constrain mean, set rotation=yes to constrain rotation.

c  If the rotation is set then a rotation angle must be set as follows:

c  rotation angle=7.0

c

c=====================================

c

c  Reference candidate (or party) name

c

c=====================================

c

c  Reflect axis=one (or two or one & two) to reflect the candidate/party map

c  through the axis. You can type either 'axis' or axes".

c  Two examples:

c    reflect axes=one & two   */ Reflect axes

c    reflect axes=1 & 2       */ Reflect axes

c  Leave the line blank up to the delimitor * for no reflections or use

c   reflect axes=             */ Reflect axes

c  The default is a blank line - no reflections

c

c=====================================

c

c   Missing data=%  for missing data flag. Example: missing data=99.

c   Any row that has one or more %'s will be not be read. That resondents scores

c    will not be used.

c

c=====================================

c

c   Bootstraps=%

c   The data is resampled % times with replacement to obtain an estimate of the

c   sampling properties of the estimates of the mean ideal point, rotation, and

c   standard deviations of the ideal points using the standard bootstrap method.

c

c=====================================

c

c  Respondent map flags=%,%,%,%,%  (1 <= 7 non-negative integer flags)

c

c  Type seven (maximum) non-negative integers separated by commas after the

c  symbol = to select the respondents who are used to estimate the political

c  map. The word 'respondent' can be omitted.

c

c  Example: map flags=1,2,10   (no comma after the last flag)

c  Separate the integers by commas with no comma after the last integer.

c

c  These flag integers are placed in the FIRST column to the right of the

c  columns of scores.

c  If a respondent's flag type in the first column equals one of the flags

c  on this line then that respondent's grades are used.

c  To omit a set of respondents whose flag is say 1 do not use 1 for any of

c  the flags.

c

c  Type the string 'use all respondents' on the line. Then all the respondent's

c  scores will be used. The program checks to see if the string 'usa all' is

c  present on a line in map.cnl above the # delimiter.

c

c  IMPORTANT: Do not type the symbol ' in the string. The symbol ' is used

c  to set off the character string from the instruction sentences above.

c  For example: use all respondents is the correct string to use all

c  repsondents.

c

c  If there is no column of respondent map flags in the data file then you

c  must type the string 'use all' to stop the program from trying to

c  find that column.

c

c  Respondent plot flags=%,%,%,%,%  (1 <= 7 non-negative integer flags)

c  Type seven (maximum) non-negative integers after the symbol = to select

c  respondents for ideal point calculations from the respondents selected by

c  either the map flags or the use all command.

c  Separate the integers by commas with no comma after the last integer.

c

c  Example: plot flags=1,2,7

c

c   These flag integers are place in the SECOND column to the right of the

c   columns of scores.

c   The selection method is the same as for the input of respondent scores

c   but the flags integers can be different.

c

c   If there is no column of ideal point flags in the data file then you must

c   type the string 'plot all' to stop the program from trying to find

c   that column.

c=====================================

c

c  # */ Candidate or party list with valence scores if available

c  Names of the m candidates/parties underneath the above title

c

c  If there are valence scores type them as v=% on the same line as the

c  candidate/party name as for example v=1 for candidate 1.

c

c  Type 'skip' or 'no' to the right of a candidate/party you do not want

c  to use. The program will rearrange the data for you.

c

c  Type 'vote' to the right of a candidate if you want to estimate the

c  proportion of the respondents who will vote for that candidate/party.

c  Type at least ONE vote in the list.

c

c  IMPORTANT: Place the symbol * after each name. The program searches for

c  the symbol * to find the end of the name, valence & vote input.

c  Also place the symbols # *  in a line after these names are read.

c  The program searches for the symbol # to find the end of the list.

c  The delimiter # ends the input of the candidate or party list.

c

c============================================

c

c  Organization of the data file structure

c

c  =======================================

c  Line 1: Name of the data

c  Line 2: n, m, format

c  n - number of respondents (or regions)

c  m - number of candidates (or parties

c  Line 3: Names of the candidates/parties over the grade columns. The order

c  of the names must match the order of the names in set in map.cnl. The

c  program checks the order on the data file against the order in the cnl file.

c

c  The data file has n rows after the three line header.

c

c  If there is respondent identification code it should be in Column 1.

c  The respondent id is read as a character string from this column if

c  the switch id=yes is set in map.cnl.

c  If there is no respondent id code column then the first column is

c  the scores for the first candidate/party.

c

c  The candidate/party grades are in m adjacent columns of the same

c  width. The grades must be right justified. If is vital that the

c  the grades columns have the same width.

c  A missing grade is coded as % where % is chosen by the user.

c   Any respondent who has at least one % is NOT USED.

c

c  RESPONDENT SELECTION FLAG COLUMNS:

c

c  Two flag columns to the right of the score columns if flags are used.

c  Note the order of these columns!

c

c  The first flag column contains integers that select the respondents whose

C  scores are to be used to obtain a map.

c  The second flag column contains integers that classify the respondents

c   whose ideal points are to be graphed.

c

c  IMPORTANT: The flag files MUST have the i format since the symbol 'i'

c  is used by the program to identify these flag columns, and the

c  scores must be read by the f (fixed) format.

c

c  Format - enclosed in ( ) for the scores data array in 50 characters or less.

c  Suppose that there are ten candidates with valence scores and a

c  respondent five character id column and two respondent flag columns

c  whose maximum value is 99. Assume that the two flag columns have a

c  one space separation. Also assume that the candidate grades have a

c  maximum value of 100 and a minimum value of 0 and that all columns

c  are separated by one space.

c

c  The format on line 2 of the data file must be of the form

c  (a5,20(1x,f3.0),2(1x,i2)). The a5 format calls for a five character

c  respondent id. The f3.0 calls for a number of at most three spaces

c  (such as 100) with nothing to the right of the decimal point or no

c  decimal point. The 1x is a one space separation between the numbers

c  on each row. The i2 format calls for integers.

c  If the grades are between 0 & 10 use (a5,20(1x,f2.0),2(1x,i2)).

c  Another example is (10(f7.1),i2,i2) if one is using % aggregate data for

c  ten candidates or parties and two adjacent integer flag columns.

c

c  A format (10(f7.1),i1,i2) for the data structure requires that each

c  row contain ten real numbers of the form -12345.6  or 12345.6 lined up

c  with one value to the right of the decimal and an 11th column with two

c  adjacent integers from 0 to 9, such as 17. It is better to use 1x

c  one space separators in order to read the data file.

c

c  If the first column contains lables for the respondents that are only

c  three characters wide then use a3 instead of a5.

c

c  If there are no candidate valence scores for ten candidates and the

c  two flag columns and grades from 0 to 100 then use (a5,10(1x,f3.0),2(1x,i2))

c  or the equivalent integer i format.

c

c  If there is only one flag column then use (a5,20(1x,i3),1x,i2)

c  for ten candidate grades and ten valence score columns.

c

c  Any book or manual on Fortran F77 or F90 provides examples of f and i formats.

c  Always use the same format for each column.

c=====================================

c  Number of candidates required to make a spatial map for dim = 2

c  Valence scores for candidates => 7    No valunce score in the data => 6

c============================================

      allocatable::can,cand,check,uc,uv,va

      real va(:)

      integer uc(:),uv(:),flgd(7),flgi(7),vi,pow,dim,pr,cns

      character*700 name

      character*80 buf

      character*50 charead,fmt

      character*30 can(:),cand(:),check(:),csub

      character*20 par

      character*2 delim,ref

      character dav

      logical exs

      real r

c********************************************

c   Open control and summary output file

      open(1,file='map.cnl',err=1,iostat=io1,status='old')

      open(3,file='Map.out',err=2,iostat=io2)

c  Start time

      call cpu_time(gs)

c-------

c  Write name of program

      write(3,'(2x,''The Cahoon-Hinich Scaling Algorithm for a Spatial A

     *nalysis of Candidate Score'')')

      ip=0

      do k=1,21

c  Clear buf

       do n=1,80

        buf(n:n)=' '

       enddo

c  Read line

       read(1,'(a80)') buf

       ia=index(buf,'#')

       if(ia>0) then

        exit

       endif

       ig=index(buf,'*')

       if(ig==0) then

        write(3,*) 'Place comment marker * after parameters on map.cnl'

        write(3,'(a80)') buf

        stop

       endif

       ig=ig-1

       name(ip+1:ip+ig)=buf(:ig)

       ip=ip+ig

      enddo

c  Data file name from control file

      par='ile name'

      delim(1:2)='= '

      fmt=charead(par,name,delim,ia,ier,3)

      if(ier>0) then

       write(3,*) 'Error in reading the file name from map.cnl'

       write(3,*) 'Error number for read = ',ier

       stop

      endif

c-------

      close(2)

c  Inquire if file exists

      inquire(file=fmt,exist=exs)

      if(.not.exs) then

       write(3,*) 'Input file ',fmt,' does not exist'

       stop

      endif

      open(2,file=fmt,err=3,iostat=io3,status='old')

      rewind(2)

c  Read & write file identity

      read(2,'(a80)') buf

      write(3,'(a80)') buf

c  Sample size in data file, no. of candidates

      read(2,*) n,m

      write(3,'(/7x,'' Number of respondents = '',i7)') n

      write(3,'(7x,'' Number of candidates in the data file = '',i3/)')

     *m

c  Euclidean distance power, dimension

      ier=0

      par='distance'

      fmt=charead(par,name,delim,ia,ier,3)

      i1=index(fmt,'stra')

      i2=index(fmt,'squar')

      if(i1>0) then

       pow=2

        elseif(i2>0) then

         pow=1

          else

           write(3,*) 'Distance power switch is wrong'

           write(3,*) 'Check spelling in ',fmt(1:30),' on map.cnl'

           stop

      endif

      par='dimension'

      fmt=charead(par,name,delim,ia,ier,3)

      i1=index(fmt,'one')

      i1n=index(fmt,'1')

      i2=index(fmt,'two')

      i2n=index(fmt,'2')

      if(i1>0 .or. i1n>0) then

       dim=1

        elseif(i2>0 .or. i2n>0) then

         dim=2

          else

           write(3,*) 'Political space dimension is not one or two'

           write(3,*) 'Check spelling in map.cnl'

           stop

      endif

c  Valence dimension, voter id

      vi=0

      par='nce dimension'

      fmt=charead(par,name,delim,ia,ier,3)

      if(ier==0) then

       i1=index(fmt,'yes')

       if(i1>0) then

        vi=1

       endif

      endif

      if(ier>0) then

       write(3,*) 'Error in read of valence dimension switch ',ier

       stop

      endif

c  Voter id

      dav='n'

      par='oter id'

      fmt=charead(par,name,delim,ia,ier,3)

      i1=index(fmt,'yes')

      if(i1>0) then

       dav='y'

      endif

c  Plot map & ideal points

      par='graph files'

      fmt=charead(par,name,delim,ia,ier,3)

      i2=index(fmt,'yes')

      if(i2>0) then

       pr=1

        else

         pr=0

      endif

      if(pow==2) then

       write(3,'(2x,''Scores are Linear in Euclidean Distance'')')

        else

         write(3,'(2x,''Scores are Linear in Squared Euclidean Distance'

     *')')

      endif

      if(dim==2) then

       write(3,'(/7x,''Political Space is Two Dimensional'')')

        else

         write(3,'(/7x,''Political Space is One Dimensional'')')

      endif

c-------

c  Check for switches for mean or rotation contraints

      cns=0

      par='ean'

      fmt=charead(par,name,delim,ia,ier,3)

      i1=index(fmt,'yes')

      i2=index(fmt,'YES')

      if(i1>0 .or. i2>0) then

c  Constrain mean to 0

       cns=1

      endif

      par='otation'

      fmt=charead(par,name,delim,ia,ier,3)

      i1=index(fmt,'yes')

      i2=index(fmt,'YES')

      if(i1>0 .or. i2>0) then

c  Constrain rotation

       cns=2

      endif

      if(i1>0 .or. i2>0) then

c  Find rotation angle

       par='angle'

       r=rnumread(par,name,delim,ia,ier,io)

      endif

c  Reference candidate or party

      par='andidate'

      csub=charead(par,name,delim,ia,ier,3)

      if(ier==1) then

       par='arty'

       csub=charead(par,name,delim,ia,ier,3)

      endif

c  Trap scatter'

      if(sc<0.) then

       write(3,*) 'Scale ',sc,' < 0'

       stop

        elseif(sc>.5) then

         write(3,*) 'Scale ',sc,' > 0.5'

         stop

      endif

c  Iterations max likelihood search

      ntry=250

c------

      num=-1

c  Flags for respondents in map

      i1=index(name,'se all')

      if(i1>0) then

c  Use all respondents with -1

       do k=1,7

        flgd(k)=-1

       enddo

      endif

      if(i1==0) then

       par='ap flags'

       fmt=charead(par,name,delim,ia,ier,3)

       if(ier/=0) then

        write(3,'(/'' ERROR reading map flag integers''/)')

        write(3,'(7x,''Syntax is: respondent map flags='')')

        stop

       endif

       i1=index(fmt(:2),' ')

       if(i1==2) then

        call number(fmt(:10),1,num,3)

        flgd(1)=num

        do k=2,7

         flgd(k)=flgd(1)

        enddo

       endif

       if(i1/=2) then

        ns=0

        i1=1

        i2=1

        do while(i2>0)

         i2=index(fmt(i1:),',')

         ic=i2-1

         call number(fmt(i1:i1+9),ic,num,3)

         i1=i1+i2

         ns=ns+1

         flgd(ns)=num

        enddo

        i2=index(fmt(i1:),' ')

        ic=i2-1

        call number(fmt(i1:i1+9),ic,num,3)

        flgd(ns)=num

        if(ns<7) then

         do k=ns+1,7

          flgd(k)=flgd(ns)

         enddo

        endif

       endif

      endif

c  Voter flags for ideal points

      i1=index(name,'lot all')

      if(i1>0) then

c  Plot all respondents with -1

       do k=1,7

        flgi(k)=-1

       enddo

      endif

      if(i1==0) then

       par='lot flags'

       fmt=charead(par,name,delim,ia,ier,3)

       i1=index(fmt(:2),' ')

       if(i1==2) then

        call number(fmt(:10),1,num,3)

        flgi(1)=num

        do k=2,7

         flgi(k)=flgi(1)

        enddo

       endif

       if(i1/=2) then

        ns=0

        i1=1

        i2=1

        do while(i2>0)

         i2=index(fmt(i1:),',')

         ic=i2-1

         call number(fmt(i1:i1+9),ic,num,3)

         i1=i1+i2

         ns=ns+1

         flgi(ns)=num

        enddo

        i2=index(fmt(i1:),' ')

        ic=i2-1

        call number(fmt(i1:i1+9),ic,num,3)

        flgi(ns)=num

        if(ns<7) then

         do k=ns+1,7

          flgi(k)=flgi(ns)

         enddo

        endif

       endif

      endif

c  Reflect axis

      ia=0

      ref='nn'

      par='eflect axes'

      fmt=charead(par,name,delim,ia,ier,3)

      if(ier==1) then

       par='eflect axis'

       fmt=charead(par,name,delim,ia,ier,3)

      endif

      if(ia>0) then

       buf(:12)=name(ia:ia+12)

       i1=index(buf(:12),'one')

       i2=index(buf(:12),'1')

       if(i1>0 .or. i2>0) then

        ref(1:1)='y'

       endif

       i1=index(buf(:12),'two')

       i2=index(buf(:12),'2')

       if(i1>0 .or. i2>0) then

        ref(2:2)='y'

       endif

      endif

      if(ref(1:1)=='y' .and. ref(2:2)=='n') then

       write(3,'(/7x,''Map is reflected through axis one'')')

      endif

      if(ref(1:1)=='n' .and. ref(2:2)=='y') then

       write(3,'(/7x,''Map is reflected through axis two'')')

      endif

      if(ref(1:1)=='y' .and. ref(2:2)=='y') then

       write(3,'(/7x,''Map is reflected through axes one & two'')')

      endif

c-------

c   Missing data flag

      par='ing data'

      mis=rnumread(par,name,delim,ia,ier,io)

      smis=float(mis)

      if(ier>0) then

       smis=numread(par,name,delim,ia,ier,io)

       if(ier>0) then

        write(3,*) 'Error in read of missing data flag'

        stop

       endif

      endif

c  Resamples

      par='straps'

      nstr=numread(par,name,delim,ia,ier,3)

      if(ier>0) then

       write(3,*) 'Error in read of bootstraps'

       stop

      endif

      if(nstr<0) then

       write(3,*) 'Error in read of bootstraps - ',nstr

       stop

        elseif(nstr==0) then

         nstr=1

      endif

c-------

c   Data format

      backspace(2)

      read(2,'(a80)') name

      ia=index(name,'(')

      ib=ia-1+index(name(ia:),')')

      ic=index(name(ib+1:),')')

      if(ic>0) then

       ib=ib+ic

      endif

      ic=index(name(ib+1:),')')

      if(ic>0) then

       ib=ib+ic

      endif

      ic=index(name(ib+1:),' ')

      if(ib>ia .and. ic>0) then

       ipar=ib-ia+1

       if(ipar>50) then

        write(3,*) 'Format is longer than 50 characters'

        pause

        stop

       endif

       fmt(1:(ipar))=name(ia:ib)

        else

         write(3,*) 'Format on second line of data file is incorrect'

         pause

         stop

      endif

c  Trap /= id & flag formats

      ia=index(fmt,'i')

c  Trap integer format for the last column of integer flags

      if(ia>0) then

       ib=index(fmt(ia:),'i1')

       ic=index(fmt(ia:),'i2')

       id=index(fmt(ia:),'i3')

       if(ib==0 .and. ic==0 .and. id==0) then

        write(3,'(/'' Use an i1, i2 or i3 integer format on line 2 of th

     *e data file'')')

        write(3,'('' to control the reading of voter integer flags''/)')

        write(3,'(a50)') fmt

        write(3,'(/'' Correct the format for the last column with the in

     *teger flags'')')

        stop

       endif

      endif

c  Trap no i in fmt if flags are set

      if(ia==0 .and. flgd(1)>=0) then

       write(3,'(/'' There is no i format for the flag column in the for

     *mat ''/,a50)') fmt

       stop

      endif

      if(dav=='y') then

       ia=index(fmt,'a')

       ib=index(fmt,'(a')

       if(ia==0) then

        write(3,'('' There is no a in the format on line 2 of the data f

     *ile to read the voter id lable''/)')

        write(3,'(a50)') fmt

        write(3,'(/'' Correct the format for the first column with the v

     *oter indentification values'')')

        stop

         elseif(ib==0) then

          write(3,'('' The a character format in the format on line 2 of

     * the data file is not of the form (a##,''/)')

          write(3,'(a50)') fmt

          write(3,'(/'' Correct the a format'')')

          stop

       endif

      endif

      if(pr==1) then

c Output ideal points to file 7

       open(4,file='Choices.graph',err=4,iostat=io4)

       open(7,file='Voters.graph',err=5,iostat=io5)

      endif

c-------

c  Allocate space

      m1=m-1

      allocate(can(0:m1),cand(0:m1),check(0:m1),uc(0:m1),uv(0:m1),va(0:m

     *1),stat=ibad)

      if(ibad/=0) then

       write(io,*) 'Unable to allocate work space for arrays in Map'

       stop

      endif

      call vot(can,cand,uc,uv,va,flgd,flgi,m,n,pow,dim,pr,csub,dav,cns,v

     *i,fmt,ref,sc,r,smis,ntry,check,nstr,3)

      write(3,*)

c+++++++++++++++++++++++++++++++++++++++++++

c  Time and date of run & run time

      call head(mt,ndy,ny,gs,rt,3)

c+++++++++++++++++++++++++++++++++++++++++++

      deallocate(can,cand,check,uc,uv,va)

      stop

c-------

  1   write(6,*) '**** Error on opening control file 1 ',io1

      stop

  2   write(6,*) '**** Error on opening output file ',io2

      stop

  3   write(3,*) '**** Error on opening data file ',io3

      stop

  4   write(3,*) '**** Error on opening candidates/parties file ',io4

      stop

  5   write(3,*) '**** Error on opening voter ideal point file ',io5

      stop

      end

c

c===========================================

c

      subroutine vot(can,cand,uc,uv,va,flgd,flgi,m,n,pow,dim,pr,csub,dav

     *,cns,vi,fmt,ref,sc,r,smis,ntry,check,nstr,io)

c********************************************

c  m - no. candidates in list, mu - no. used

c============================================

      allocatable::a,c,cr,cov,d,dbar,e,er,f,id,s,typ,v,vs,vd,w,x

      real x(:),c(:),cr(:),va(m),vs(:)

      real*8 a(:),cov(:),d(:),dbar(:),v(:),er(:),f(:),e(:),s(:),w(:)

      integer id(:),typ(:),uc(m),uv(m),flgd(7),flgi(7),vi,pow,dim,pr,psu

     *b,cns

      character*700 name

      character*80 buf

      character*50 fmt,fmo

      character*30 csub

      character*30 can(m),cand(m),check(m)

      character*20 par

      character*5 vd(:)

      character*2 delim,ref

      character dav

      intent(in) flgd,flgi,m,n,pow,dim,pr,csub,cns,vi,fmt,ref,sc,r,smis,

     *ntry,nstr,io

c********************************************

c  Read candidate names

      mu=0

      psub=0

      nsum=0

      vmax=0.

      do k=1,m

       read(1,'(80a)') name

c  Check for blanks

       ia=index(name(1:1),' ')

       ib=index(name(2:2),' ')

       ic=index(name(3:3),' ')

       if(ia>0 .and. ib>0 .and. ic>0) then

        write(io,'('' There are 3 blank columns in the candidate/party n

     *ame in Map.cnl'')')

        write(io,'(a77)') name(:77)

        stop

       endif

       ia=index(name,'#')

       if(ia>0) then

        if(k-1/=m) then

         write(io,'(/''ERROR''/)')

         write(io,'('' The number of candidates '',i2,'' set in the data

     * file is not the same as the'',/,'' number of candidates in map.cn

     *l'')') m

         stop

        endif

        exit

       endif

       ig=index(name,'*')

       if(ig==0) then

        write(io,*) 'Place comment marker * after candidate/party names

     *and vote switches in the list'

        write(io,'(a77)') name(:77)

        stop

       endif

       va(k)=0.

       uc(k)=0

       ia=index(name,'=')

       if(ia>0) then

        cand(k)=name(:ia-3)

         else

c  Find blank after character name

          ib=index(name,' ')

          cand(k)=name(:ib)

       endif

       i1=index(name(4:),' no ')

       i2=index(name(4:),' skip ')

       if((i1>0.or.i2>0).and.cand(k)==csub) then

        write(io,'(2x,''ERROR - you are skipping the reference candidate

     * '',i2,1x,a30)') k,csub

        stop

       endif

       if(i1==0 .and. i2==0) then

c  Candidate/party k is mapped

        uc(k)=1

        mu=mu+1

        can(mu)=cand(k)

        if(can(mu)==csub) then

         psub=mu

        endif

c  Find vote count flag

        ia=index(name,' vote ')

        uv(mu)=0

        if(ia>0) then

c  Compute candidate vote count

         uv(mu)=1

         nsum=nsum+1

        endif

        delim='='

        if(vi==1) then

c  Find valence value if v= is present

         par='v'

         ve=rnumread(par,name,delim,ia,ier,io)

         va(mu)=ve

         if(ier>0) then

          write(io,*) 'Error in read of valence ',k,' ',par

          stop

         endif

         if(ve>vmax) then

          vmax=ve

         endif

        endif

       endif

      enddo

c  End of candidate loop

      do k=1,mu

c  Highest valence is 0

       va(k)=vmax-va(k)

      enddo

c-------

c  Trap psub

      if(psub==0) then

       write(io,*) 'Check spelling of the name of the reference candidat

     *e'

       write(io,*) csub

       stop

      endif

c  Trap nsum

      if(nsum==0) then

       write(io,'(/'' No candidate or party vote is being calculated'')'

     *)

       write(io,'('' Type vote to the right of a candidate or party to g

     *et a test vote'')')

       stop

      endif

c  Read candidate name line

      read(2,'(a700)') name

      call lineread(name,check,m ,io)

      do kc=1,m

       if(cand(kc)/=check(kc)) then

        write(io,'(/''ERROR''/)')

        write(io,'('' Candidate-Party Name '',a30,/'' is not the same as

     * the name over column '',a30)') cand(kc),check(kc)

        write(io,'(/7x,''Verify that the number of candidates '',i2,'' s

     *et in the data file'',/,'' is the same as the number of candidates

     * in map.cnl''/)') m

        write(io,'(2x,''Verify that there is at least one space between

     * the candidate names in the 3rd line of the data file'')')

        stop

       endif

      enddo

c-------

      if(vi==1) then

       write(io,'(/7x,''Candidate Valence Scores Used'')')

      endif

      write(io,'(/7x,a20,'' scores are subtracted from the others'')') c

     *an(psub)

      if(cns==1 .or. cns==3) then

       write(io,'(7x,''This candidate position is set at the origin'')')

      endif

c-------

      do k=1,50

       fmo(k:k)=' '

      enddo

      if(mu==m .and. dav=='n') then

c  No respondent id in first column & all candidates are mapped

       ia=index(fmt,'a')

       if(ia>0) then

c  Replace a# by #x, to remove id character read

        ib=index(fmt(ia+1:),',')

        fmo(1:1)='('

        fmo(2:ib)=fmt(ia+1:ia+ib-1)

        nt=ib+1

        fmo(nt:nt)='x'

        fmo(nt+1:)=fmt(ia+ib:)

         else

c  No a in fmt, shift to fmo

          fmo(:50)=fmt(:50)

       endif

      endif

c-------

c  Respondent id in first column & all candidates are mapped

      if(mu==m .and. dav=='y') then

c  Shift to fmo with voter id format

       fmo(:50)=fmt(:50)

      endif

c  No respondent id in first column & some candidates are skipped

      if(mu<m .and. dav=='n') then

       ia=index(fmt,'a')

       if(ia>0) then

c  Replace a# by #x, to remove id character read

        ib=index(fmt(ia+1:),',')

        fmo(1:1)='('

        fmo(2:ib)=fmt(ia+1:ia+ib-1)

        nt=ib+1

        fmo(nt:nt)='x'

        fmo(nt+1:)=fmt(ia+ib:)

         else

c  No a in fmt, shift to fmo

          fmo(:50)=fmt(:50)

       endif

      endif

c  Respondent id in first column & some candidates are skipped

      if(mu<m .and. dav=='y') then

       fmo(:50)=fmt(:50)

      endif

c-------

      if(dim==2) then

c  Required no. of candidates

       if(vi==1) then

        if(mu<8) then

         write(io,*) 'Number of candidates ',mu,' < 8. Use more candidat

     *es'

         stop

        endif

       endif

       if(vi==0) then

        if(mu<6) then

         write(io,*) 'Number of candidates ',mu,' < 7. Increase no. of c

     *andidates used'

         pause

         stop

        endif

       endif

      endif

      if(cns==2) then

         write(io,'(/7x,''Rotation set to angle '',f7.1)') r

      endif

      write(io,'(/7x,''Number of candidates used = '',i3)') mu

      write(io,'(/7x,''Data format used = '',a50)') fmo

      write(io,'(/7x,''Missing data flag = '',f6.1)') smis

      write(io,'(/7x,''Ideal Point Scatter Parameter = '',f7.4)') sc

c-------

      m1=m-1

      mn=mu-1

      nn=n-1

      nmu=n*mu-1

      mu2=mu*mu-1

      mn2=mn*mn-1

      allocate(x(0:nmu),id(0:nn),typ(0:nn),cov(0:mn2),dbar(0:mn),a(0:mu2

     *),d(0:mn2),v(0:mn2),er(0:mn),f(0:2*m),e(0:mn),w(0:mu2),c(0:nmu),cr

     *(0:2*mu-1),s(0:mn2),vs(0:m1),vd(0:nn),stat=ibad)

      if(ibad/=0) then

       write(io,*) 'Unable to allocate work space for arrays in top'

       stop

      endif

      call top(can,cand,a,c,cov,cr,d,dbar,e,er,f,id,s,typ,uc,uv,v,va,vs,

     *vd,w,flgd,flgi,m,mu,mn,n,pow,dim,pr,psub,cns,dav,vi,fmo,ref,sc,r,s

     *mis,nstr,ntry,io)

c------

      deallocate(a,c,cr,cov,d,dbar,e,er,f,id,typ,v,vs,vd,w,x)

      return

      end

c

c===========================================

c

      subroutine top(can,cand,a,c,cov,cr,d,dbar,e,er,f,id,s,typ,uc,uv,v,

     *va,vs,vd,w,flgd,flgi,m,mu,mn,n,pow,dim,pr,psub,cns,dav,vi,fmo,ref,

     *sc,r,smis,nstr,ntry,io)

c********************************************

c  uc - candidate selection 0/1, uv - vote count 0/1, va - valences

c============================================

      allocatable::iw,us

      real x(n,mu),c(n,mu),cr(mu,2),va(m),vs(m),y(7),bmn(4),br(2),bs1(2)

     *,bs2(2),bv(2)

      real*8 a(mu*mu),cov(mn,mn),d(mn,mn),dbar(mu),e(mu),er(mu),f(mn),s(

     *mn*mn),v(mn,mn),w(mn*mn),z1

      integer id(n),iw(:),uc(m),us(:),uv(m),typ(n),flgd(7),flgi(7),vi,po

     *w,dim,pr,psub,cns,t

      character*50 fmo

      character*30 can(m),cand(m)

      character*5 vd(n)

      character*2 ref

      character dav

      intent(in) can,flgd,flgi,uc,uv,va,m,mu,mn,n,pow,dim,pr,cns,vi,fmo,

     *ref,sc,r,smis,nstr,ntry,io

      intent(out) cand,vs

c********************************************

      allocate(iw(0:mn),us(0:mn),stat=ibad)

      if(ibad/=0) then

       write(io,*) 'Unable to allocate work space for iw in run'

       stop

      endif

c-------

c  Rearrange data x =>c & valence va =>vs

      call dat(x,typ,dbar,can,cand,va,vs,vd,uc,us,uv,e,c,a,m,mu,n,nu,pow

     *,psub,smax,smin,fmo,flgd,flgi,dav,smis,vi,io)

      write(3,'(/7x,''Maximum & Minimum Values for the Scores: '',f7.1,2

     *x,f7.1)') smax,smin

c  Trap nu

      if(nu<1) then

       write(io,'(/7x,''Error'')')

       write(io,*) 'No respondents were selected for the flags set in ma

     *p.cnl'

       stop

      endif

      if(flgd(1)>=0) then

       write(io,'(/7x,i5,'' respondents for estimating the map were sele

     *cted using the flags:''/)') nu

       write(io,'(7('' Map('',i2'')= '',i2))') (k,flgd(k),k=1,7)

        else

         write(io,'(/7x,'' All '',i7,'' respondents were used to map the

     * space'')') nu

      endif

c  Covariance matrix of adjusted scores

      call covm(c,cov,dbar,n,nu,mn,io)

c-------

      call maxle(cov,nu,mn,dim,a,d,e,er,f,v,w,sge,z1,ntry,1,ier,io)

      er(mu)=dble(sge*sge)

c  Error standard deviations

      write(io,'(/7x,''Specific Error Standard Deviations''/)')

      write(io,'(2x,a20,2x,f15.4)') (cand(k),dsqrt(er(k)),k=1,mu)

c  Eigenvalues of (2/n)*I**-1

      scale=2./float(n)

      call eigen(a,mn,mn,f,e,s,jev)

      semax=sqrt(scale/f(mn))

      semin=sqrt(scale/f(1))

      write(io,'(2x,72(''=''))')

      write(io,'(/7x,''Max Standard Errors of Estimates from the Eigenva

     *lues'',/,'' of the Covariance Matrix of the Maximum Likelihood Est

     *imates = '',e12.3)') semax

      write(io,'(/7x,''Min Standard Errors of Estimates from the Eigenva

     *lues'',/,'' of the Covariance Matrix of the Maximum Likelihood Est

     *imates = '',e12.3)') semin

      write(io,'(/7x,''Ratio Max/Min = '',f12.3)') semax/semin

      write(io,'(2x,72(''=''))')

      write(io,'(7x,''OLS Estimates of Remaining Parameters''/)')

      if(dim==2) then

       call strap(a,c,cov,d,dbar,e,er,f,s,v,w,mu,mn,n,nu,ntry,vi,vs,cns,

     *dim,bmn,br,bs1,bs2,bv,nstr,io)

       ro=br(1)

       call v2(mu,v,cr,cns,r,ro,bs1(1),bs2(1),io)

        elseif(dim==1) then

         call v1(m,v,cr,s,d,dbar,w,er,f,vi,vs,wv,cns,b1,ier,io)

      endif

c-------

      write(io,'(2x,65(''=''))')

c  Write parameters and candidate positions

      if(dim==2) then

       write(io,'(/2x,''Ideal Point Standard Deviations = '',2(f14.2,1x)

     *)') bs1(1),bs2(1)

       if(vi>0) then

        write(io,'(2x,''Candidate Characteristic Value Weight = '',f14.2

     *)') bv(1)

       endif

       pi=3.141592

       if(cns==0.or.cns==1) then

        write(io,'(2x,''Rotation = '',f7.3/)') 180.*ro/pi

         elseif(cns==2) then

          write(io,'(2x,''Rotation = '',f7.3/)') r

           else

            write(io,'(/7x,''ERROR - constraint is wrong'')')

            write(io,'(/7x,''cnl = '',i1)') cns

            stop

       endif

       write(io,'(17x,''Positions of the Candidates''/)')

       do i=1,mu

        if(ref(1:1)=='y') then

         cr1=-1.*cr(i,1)

          else

           cr1=cr(i,1)

        endif

        if(ref(2:2)=='y') then

         cr2=-1.*cr(i,2)

          else

           cr2=cr(i,2)

        endif

        write(io,'(4x,a20,4x,f12.3,2x,f12.3)') cand(i),cr1,cr2

       enddo

       write(io,'(8x,''Mean Ideal Point'',4x,f12.3,2x,f12.3)') bmn(1),bm

     *n(2)

      endif

c------

      if(dim==2.and.nstr>1) then

       write(io,'(2x,65(''*''))')

       write(io,'(/15x,''Bootstrapped Parameter Means''/)')

       write(io,'(12x,''Rotation = '',f12.4)') 180./pi*br(1)

       write(io,'(7x,''Mean Ideal Point Coordinates = '',2(f12.4,2x))')

     *bmn(1),bmn(2)

       write(io,'(4x,''Ideal Point Standard Deviations = '',2(f12.4,2x))

     *') bs1(1),bs2(1)

       if(vi>0) then

        write(io,'(7x,''Valence Weight = '',f12.4)') bv(1)

       endif

       write(io,'(/15x,''Bootstrapped Standard Errors''/)')

       write(io,'(12x,''SE of Rotation = '',f12.4)') 180./pi*br(2)

       write(io,'(7x,''SE of Mean Ideal Point Coordinates = '',2(f12.4,2

     *x))') bmn(3),bmn(4)

       write(io,'(4x,''SE of Ideal Point Standard Deviations = '',2(f12.

     *4,2x))') bs1(2),bs2(2)

       if(vi>0) then

        write(io,'(7x,''SE of Valence Weight = '',f12.4)') bv(2)

       endif

      endif

      if(dim==1) then

c  Write parameters and candidate positions for one dimension

       write(io,'(/2x,''Ideal Point Standard Deviation = '',f14.2)')b1

       if(vi>0) then

        write(io,'(2x,''Candidate Characteristic Value Weight = '',f14.2

     *)') wv

       endif

       write(io,'(/17x,''Positions of the Candidates''/)')

       do k=1,m

        write(io,'(4x,a20,4x,f12.3)') cand(k),cr(k,1)

       enddo

      endif

c-------

c  Plot ideal points

      if(dim==2) then

c  Cutoff for lowest score

       chk=abs(smin-smax)

       if(pow==2) then

        chk=chk*chk

       endif

c  Output candidates/parties to file 4

       do k=1,mu

        if(ref(1:1)=='y') then

         cr(k,1)=-1.*cr(k,1)

          else

           cr(k,1)=cr(k,1)

        endif

        if(ref(2:2)=='y') then

         cr(k,2)=-1.*cr(k,2)

          else

           cr(k,2)=cr(k,2)

        endif

        if(pr==1) then

         write(4,'(a20,1x,2(f10.2,1x))') cand(k),cr(k,1),cr(k,2)

        endif

       enddo

       call ideal2(n,nu,mu,c,id,typ,uv,cr,dbar,va,bv(1),a,e,iw,y,flgi,ch

     *k,vd,dav,pow,ref,sc,ni,pr,io)

       if(pr==1) then

        write(4,'(''Ideal Point Mean'',5x,2(f10.2,1x))') bmn(1),bmn(2)

        write(4,'(''xbar'',17x,2(f10.2,1x))') y(1),y(2)

       endif

c  Voting % estimates

       write(io,'(7x,''% of Ideal Points that are Closest to Each Candid

     *ate or Party''/)')

       do k=1,mu

        if(uv(k)==1) then

c  Vote %

         write(io,'(a20,1x,f5.1)') cand(k),100.*sngl(a(k))/float(ni)

        endif

       enddo

      endif

      if(dim==1) then

       call ideal1(n,nu,mu,c,cr,typ,iw,sc,flgi,ni,io)

       if(pr==1) then

        do k=1,ni

         ncode=nint(x(k,2))

         write(7,'(f10.2,1x,i3)') x(k,1),ncode

        enddo

        do k=1,mu

         write(4,'(a20,1x,f10.2)') cand(k),cr(k,1)

        enddo

       endif

      endif

c-------

      deallocate(iw,us)

      return

      end

c

c===========================================

c

      subroutine dat(x,typ,dbar,can,cand,va,vs,vd,uc,us,uv,e,c,a,mo,m,n,

     *nu,pow,psub,smax,smin,fmo,flgd,flgi,dav,smis,vi,io)

c********************************************

c  Read data & rearrange columns & find max & min of data                  .

c  nu - number of respondents used, ncol - flag columns (0,1,2)

c  flgd - flags for respondents used, flgi - ideal point flags

c  Data from all respondents if flgd(1)=-1, map all ideal points if flgi(1)=-1

c  m - candidates used, smis - missing data flag, mo - all candidates

c  uv - vote count flag, vd - character respondent identification

c  c - respondent used scores, vs - rearranged valence scores

c  typ - flags for plotting

c============================================

      allocatable::w

      integer typ(n),uc(mo),us(m),uv(m),psub,pow,flgd(7),flgi(7),use,vi

      real x(n,m),va(mo),vs(mo),c(n,m),w(:),smax,smin

      real*8 e(m),a(2*m),dbar(m),rn

      character*50 fmo

      character*30 can(mo),cand(mo),buf

      character*5 vd(n),rvd

      character dav

      intent(in) can,mo,m,n,pow,psub,fmo,flgd,flgi,dav,smis,uc,uv,va,vi,

     *io

      intent(out) c,x,dbar,nu,smax,smin,typ,us,vd,vs

c********************************************

      allocate(w(0:mo-1),stat=ibad)

      if(ibad/=0) then

       write(io,*) 'Unable to allocate work space for w in dat'

       stop

      endif

c  Initialize

      smax=0.

      smin=100.

      nu=0

      keep=0

      ncol=0

      if(flgd(1)/=-1.or.flgi(1)/=-1) then

c  Check for one or two flag columns

       ia=index(fmo,'i')

       if(ia>0) then

        ncol=2

        ib=index(fmo(ia+1:),'i')

        ic=index(fmo(ia+1:),'))')

c  Trap one i & )) in format

        if(ib==0 .and. ic>0) then

         write(io,'(/'' Error''/)')

         write(io,'('' Remove the double )) in the format '',a50)') fmo

         write(io,'('' There is only one letter i in the format'')')

         stop

        endif

        if(ib==0 .and. ic==0) then

         ncol=1

        endif

       endif

      endif

c  Output scores for first respondent to check against data file

      write(io,'(/2x,67(''*''))')

      write(io,'(2x,''Data for First Respondent to Check Data Read'')')

c  Voter loop

      do i=1,n

       if(ncol==0) then

        nti=-1

        nxi=-1

c  Input x(i,j) for all i =>c(n,j)

        if(dav=='y') then

c  Read respondent id

         read(2,fmo) rvd,(w(j),j=0,mo-1)

         if(i==1) then

          write(io,*) rvd,(w(j),j=0,mo-1)

         endif

          else

c  Input x(i,j) for all i with no id

           read(2,fmo) (w(j),j=0,mo-1)

           if(i==1) then

            write(io,*) (w(j),j=0,mo-1)

           endif

        endif

       endif

       if(ncol==2) then

c  Input x(i,j) for all i & two flag columns (ideal select, map type)

        if(dav=='y') then

         read(2,fmo) rvd,(w(j),j=0,mo-1),nti,nxi

         if(nti==0) then

          write(io,'(/'' ERROR - the flag in the first column is 0'')')

          write(io,'('' Check to see if there is only one flag column in

     * the data set'',/'' but the format '',a50)') fmo

          write(io,'('' calls for two'')')

          stop

         endif

         if(i==1) then

          write(io,*) rvd,(w(j),j=0,mo-1),nti,nxi

         endif

         if(flgd(1)==-1) then

c  Skip respondent flag

          nti=-1

         endif

         if(flgi(1)==-1) then

c  Skip ideal select flag

          nxi=-1

         endif

        endif

        if(dav/='y') then

         read(2,fmo) (w(j),j=0,mo-1),nti,nxi

         if(i==1) then

          write(io,*) (w(j),j=0,mo-1),nti,nxi

         endif

         if(flgd(1)==-1) then

          nti=-1

         endif

         if(flgi(1)==-1) then

          nxi=-1

         endif

        endif

       endif

c-------

c  Trap one column flags

       if(ncol==1.and.(flgd(1)/=-1.and.flgi(1)/=-1)) then

        write(io,'(/''ERROR''/)')

        write(io,'('' Both map & plot flags are set but there is only on

     *e flag column in the data file'')')

        stop

       endif

       if(ncol==1.and.(flgd(1)/=-1.and.flgi(1)==-1)) then

c  Plot all

        nxi=-1

c  Input x(i,j) for all i & one flag columns (map type)

        if(dav=='y') then

         read(2,fmo) rvd,(w(j),j=0,mo-1),nti

         if(i==1) then

          write(io,*) rvd,(w(j),j=0,mo-1),nti

         endif

        endif

        if(dav/='y') then

         read(2,fmo) (w(j),j=0,mo-1),nti

         if(i==1) then

          write(io,*) (w(j),j=0,mo-1),nti

         endif

        endif

       endif

       if(ncol==1.and.(flgd(1)==-1.and.flgi(1)/=-1)) then

c  Map all

        nti=-1

c  Input x(i,j) for all i & one flag columns (plot type)

        if(dav=='y') then

         read(2,fmo) rvd,(w(j),j=0,mo-1),nxi

         if(i==1) then

          write(io,*) rvd,(w(j),j=0,mo-1),nxi

         endif

        endif

        if(dav/='y') then

         read(2,fmo) (w(j),j=0,mo-1),nxi

         if(i==1) then

          write(io,*) (w(j),j=0,mo-1),nxi

         endif

        endif

       endif

c-------

c  Max & min

       mis=0

       do k=1,mo

        sv=w(k-1)

c  Check for missing data flag

        if(sv/=smis.and.uc(k)==1) then

c  No missing data flag

         if(sv<smin) then

          smin=sv

         endif

         if(sv>smax) then

          smax=sv

         endif

          elseif(sv==smis.and.uc(k)==1) then

c  Missing data flag

           mis=mis+1

        endif

       enddo

       if(mis==0) then

        keep=keep+1

c  No missing data in row

        use=0

        do k=1,mo

         if(uc(k)==1) then

          use=use+1

          x(i,use)=w(k-1)

         endif

        enddo

c  Trap nuse

        if(use/=m) then

         write(io,*) 'Number used  ',use, '/= candidate used ',m

         stop

        endif

c--------

c  Select voters using the flgd flages

        if(nti==flgd(1).or.nti==flgd(2).or.nti==flgd(3).or.nti==flgd(4).

     *or.nti==flgd(5).or.nti==flgd(6).or.nti==flgd(7).or.nti==-1) then

         nu=nu+1

         do j=1,m

c  Move x =>c

          c(nu,j)=x(i,j)

         enddo

c  Save respondent type for ideal point

         typ(nu)=nxi

         if(dav=='y') then

c  Save respondent ident vd

          vd(nu)=rvd

         endif

        endif

         else

c  Fill out missing data row

          do j=1,m

           x(i,j)=-1.

          enddo

       endif

c  End of no missing data loop

      enddo

c  End of voter loop

      write(io,'(2x,67(''=''))')

c-------

c  Trap range

      if(smax<=0) then

       write(io,*)

       write(io,*) 'Error'

       write(io,*) 'max score ',smax,' not positive'

       stop

        elseif(smin>=smax) then

         write(io,*)

         write(io,*) 'Error'

         write(io,*) 'smin ',smin,' > smax ',smax

         stop

      endif

c  Trap nu

      if(nu==0) then

       write(io,'(/7x,''Error'')')

       write(io,*) 'Respondents selected by the first flag = 0'

       stop

      endif

c-------

c  Voter loop for columns shifts

      if(psub/=m) then

c  Rearrange columns of x =>c, candidates =>cand, valence =>vs, uv =>us

       nb=psub-1

       nc=psub+1

       buf=can(psub)

       uh=uv(psub)

       vh=va(psub)

       do i=1,n

        if(i<=nu) then

         cst=c(i,psub)

        endif

c  Rearrange nc-m columns

        do j=nc,m

         if(i<=nu) then

          c(i,j-1)=c(i,j)

         endif

         if(i==n) then

          cand(j-1)=can(j)

          vs(j-1)=va(j)

          us(j-1)=uv(j)

          if(psub>1) then

c  Shift can => cand & va =>vs

           do kc=1,nb

            cand(kc)=can(kc)

            us(kc)=uv(kc)

            vs(kc)=va(kc)

           enddo

          endif

         endif

        enddo

        if(i<=nu) then

         c(i,m)=cst

        endif

       enddo

c  End of column shift loop

       cand(m)=buf

       us(m)=uh

       vs(m)=vh

      endif

      if(psub==m) then

       do j=1,m

        cand(j)=can(j)

        vs(j)=va(j)

       enddo

      endif

c-------

c  Second voter loop

      m2=2*m

      do k=1,m

       e(k)=0.

       a(k)=0.d0

       a(m+k)=0.d0

      enddo

      do i=1,n

       if(i<=nu) then

c  Accumulate first two moments of c

        a(m)=a(m)+dble(c(i,m))

        a(m2)=a(m2)+dble(c(i,m)*c(i,m))

       endif

       do j=1,m-1

        if(i<=nu) then

c  Check smax

         if(c(i,j)>smax) then

          write(io,'(/''ERROR''/)')

          write(io,*)'x(i,j) ',c(i,j),' > max score ',smax,' i & j ',i,j

          write(io,*) 'This observation exceeds the maximum value'

          stop

         endif

         a(j)=a(j)+dble(c(i,j))

         a(m+j)=a(m+j)+dble(c(i,j)*c(i,j))

        endif

       enddo

c-------

c  Transform scores so that 0 is the best score for a candidate

c  Square if pow=2 for column m

       if(pow==1 .and. i<=nu) then

        c(i,m)=smax-c(i,m)

        e(m)=e(m)+dble(c(i,m))

       endif

c  Square if pow=2 & subtract last column

       if(pow==2) then

        if(i<=nu) then

         c(i,m)=smax-c(i,m)

         c(i,m)=c(i,m)*c(i,m)

         e(m)=e(m)+dble(c(i,m))

        endif

       endif

       do j=1,m-1

        if(pow==1 .and. i<=nu) then

         c(i,j)=smax-c(i,j)

         e(j)=e(j)+dble(c(i,j))

        endif

        if(pow==2) then

         if(i<=nu) then

          c(i,j)=smax-c(i,j)

          c(i,j)=c(i,j)*c(i,j)

          e(j)=e(j)+dble(c(i,j))

         endif

        endif

c  Subtract last squared column from each column

        if(i<=nu) then

         c(i,j)=c(i,j)-c(i,m)

        endif

       enddo

c  End of second voter loop

      enddo

c  Means of smax-score

      rn=dfloat(nu)

      dbar(m)=e(m)/rn

      do k=1,m-1

       dbar(k)=sngl(e(k)/rn)-dbar(m)

      enddo

      write(io,'(/7x,i6,'' Respondents with No Missing Scores'')') keep

      if(vi==0) then

       write(io,'(/17x,''Means & Standard Deviations of Scores''/)')

        elseif(vi==1) then

         write(io,'(/17x,''Means, Standard Deviations of Scores & Valenc

     *e Scores''/)')

      endif

      do j=1,m

       a(j)=a(j)/rn

       a(m+j)=a(m+j)/float(nu-1)-(rn/float(nu-1))*a(j)*a(j)

       if(vi==0) then

        write(io,'(1x,a20,2(3x,f12.2))') cand(j),a(j),sqrt(a(m+j))

         elseif(vi==1) then

         write(io,'(1x,a20,3(3x,f12.2))')cand(j),a(j),sqrt(a(m+j)),va(j)

       endif

      enddo

c-------

      deallocate(w)

      return

      end

c

c===========================================

c

      subroutine covm(x,cov,dbar,n,nu,mn,io)

c********************************************

c  Covariance matrix of adjusted scores used

c============================================

      real x(n,mn)

      real*8 cov(mn,mn),dbar(mn),ck1,ck2

      intent(in) x,dbar,mn,n

      intent(out) cov

c********************************************

c  Initialize

      do k1=1,mn

       do k2=1,mn

        cov(k1,k2)=0.d0

       enddo

      enddo

c  Covariance matrix

      do i=1,nu

       do k1=1,mn

        ck1=dble(x(i,k1))-dbar(k1)

        do k2=k1,mn

         ck2=dble(x(i,k2))-dbar(k2)

         cov(k1,k2)=cov(k1,k2)+ck1*ck2

        enddo

       enddo

      enddo

      do k1=1,mn

       do k2=k1,mn

        cov(k1,k2)=cov(k1,k2)/dfloat(nu)

        cov(k2,k1)=cov(k1,k2)

       enddo

      enddo

c-------

      return

      end

c

c===========================================

c

      subroutine maxle(cov,n,m,dim,a,d,e,er,f,v,w,sge,z1,ntry,sw,ier,io)

c*************************************************

c   m - Number of targets-reference (mn in call), sw=1 output eigenvalues

c   dim - Common factors,  cov = vv'+sg0*11'+er, sw=2 output all writes

c   v - matrix of factor loadings, er - Specific variances, z1 - min of -ln(L)

c   fmin - min of likelihood objective function, sge - sg0 sig estimate

c   2nd derivative matrix of S=>d, Lawley's approximation =>a

c=================================================

      implicit real*8(a-h,o-z)

      allocatable::c

      real*8 a(m*m),c(:,:),cov(m,m),d(m,m),v(m,m),w(m,m),e(m),er(m),f(m)

      real sge

      integer dim,sw,ww,t

      intent(in) cov,n,m,ntry,sw,io

      intent(out) a,d,er,e,ier,sge,v,z1

c*************************************************

      allocate(c(0:m-1,0:m-1),stat=ibad)

      if(ibad/=0) then

       write(io,*) 'Unable to allocate c in maxle'

       stop

      endif

      ier=0

c  Find minimum fac for mle of sg0

      t=0

      fmin=1.d7

      sgs=1.d0

      dowhile(t<ntry.and.ier==0)

       t=t+1

       str=1.d-3+dfloat(t)/dfloat(ntry+1)

c  cov-str*1*1' => w

       do k1=1,m

         do k2=1,k1

          cd=cov(k1,k2)-dble(str*str)

          w(k1,k2)=cd

          w(k2,k1)=cd

         enddo

        enddo

c  Factor w - no output if sw=0

       call cfac2(w,v,n,m,dim,a,d,er,e,f,c,z1,0,ier,io)

       if(ier==0) then

        if(z1<fmin) then

         fmin=z1

         sgs=str

        endif

       endif

      enddo

c  End of mle loop

      do k1=1,m

       do k2=1,k1

         cd=cov(k1,k2)-sgs*sgs

         w(k1,k2)=cd

         w(k2,k1)=cd

       enddo

      enddo

      sge=sngl(sgs)

      call cfac2(w,v,n,m,dim,a,d,er,e,f,c,z1,sw,ier,io)

      if(ier==1) then

       write(io,*) 'Cov is not positive definite in maxle'

      endif

c  Divide by 2

      do k=1,m

       v(k,1)=v(k,1)/2.d0

       v(k,2)=v(k,2)/2.d0

      enddo

c-------

      deallocate(c)

      return

      end

c

c===========================================

c

      subroutine cfac2(cov,v,nr,n,nvec,a,d,er,e,f,s,z1,sw,ier,io)

c*************************************************

c  Minimize f = sum(j=nvec+1,n)(e(j)-log(e(j)-1)

c  Input: n - Dimension of sample covariance matrix

c  nr - Number of respondents in sample used, sw=1 - output eigenvalues

c  nvec - Number of common factors, cov - Sample covariance matrix

c  Ouput: v - n x nvec matrix of factor loadings, z1 - min of -ln(L)

c  er - Specific variances, e - Eigenvalues of (er**-1/2)*cov*(er**-1/2)

c  2nd derivative matrix of S=>d, Lawley's approximation =>a, sw=2 - all writes

c=================================================

      implicit real*8(a-h,o-z)

      allocatable::indx,w

      real*8 a(n,n),cov(n,n),d(n,n),s(n,n),v(n,n),er(n),f(n),e(n),w(:),h

      real dd

      integer indx(:),sw

      intent(in) cov,n,nr,nvec,sw,io

      intent(out) a,d,er,e,v,z1,ier

c*************************************************

      allocate(indx(0:n-1),w(0:n-1),stat=ibad)

      if(ibad/=0) then

       write(io,*) 'Unable to allocate indx in cfac2'

       stop

      endif

c  Max for loop

      mxt=100

c  Initialize z1 & con2

      z1=1.0d10

      con2=cov(1,1)*7.d-2

      kp1=nvec+1

      ier=0

c-------

c  Cov => a & d

      do i=1,n

       f(i)=0.

       er(i)=0.d0

       do j=i,n

        d(i,j)=cov(i,j)

       enddo

      enddo

      inv=0

c  w=d**-1*f=0

      call dls(d,w,f,e,n,inv,ier,io)

c  Test if cov is positive definite

      if(ier==1.and.sw==2) then

       write(io,*) 'Cov is not positive definite in cfac2'

       return

      endif

c  Iniitialize er

      do i=1,n

c  er=(1-nvec)/(2*n*diag(cov))

       er(i)=(1.d0-dfloat(nvec)/(2.d0*dfloat(n)))/d(i,i)

      enddo

c  Eigenvalues of (er**-1/2)*cov*(er**-1/2) => e, eigenvectors => v

      call phi(cov,a,v,e,er,w,n,nvec,io)

c  1st derivative of objective S => f

      call fd(cov,v,e,er,n,nvec,f)

c  w=a**-1*f

      call dls(a,w,f,e,n,0,ier,3)

      if(ier==1.and.sw==2) then

       write(io,*) 'a is not positive definite'

       ier=2

       return

      endif

c  er-w => er

      do i=1,n

       er(i)=er(i)-w(i-1)

       if(er(i)<1.d-10) then

        er(i)=1.d-10

       endif

      enddo

c-------

      ni=0

      z2=0.

c  Start of MLE loop

      dowhile (z2<0.99999 .and. ni<=mxt)

       ni=ni+1

       do i=1,n

        if(er(i)<1.d-10) then

         er(i)=1.d-10

        endif

       enddo

c  Eigenvalues of (er**-1/2)*cov*(er**-1/2) => e, eigenvectors => v

       call phi(cov,a,v,e,er,w,n,nvec,io)

       if (z1<=0.d0) then

        z1=1.d-2

       endif

c  1st derivatives of the log likelihood S => f

       call fd(cov,v,e,er,n,nvec,f)

c  2nd derivative matrix of the log likelihood S => d

       call sd(v,e,er,n,nvec,f,a,d)

c  Objective function Sum(e-log(e)-1)

       z=0.d0

       do k=kp1,n

        z=z+e(k)-dlog(e(k))-1.d0

       enddo

       z2=z1

       z1=z

c  f=>e & d=>s

       do k=1,n

        e(k)=f(k)

        s(k,k)=d(k,k)

        do j=k+1,n

         h=d(j,k)

         s(j,k)=h

         s(k,j)=h

        enddo

       enddo

c  w=d**-1*f

       call dludcmp(s,n,n,indx,dd,ier,io)

       call dlubksb(s,n,n,indx,w)

       nchk=0

       do j=1,n

        if(dabs(w(j-1))>con2) then

         nchk=nchk+1

        endif

        if(nchk==1) then

         do k=1,n

          w(k-1)=e(k)

         enddo

        endif

       enddo

       do k=1,n

        er(k)=er(k)-w(k-1)

       enddo

       z2=z1/z2

      enddo

c  End of MLE loop

c-------

      if(ni==mxt) then

       write(io,'(''Loop terminated at max '',i4,'' iterations'')') mxt

       ier=3

       return

      endif

c  Trap negative er's

      do k=1,n

       if(er(k)<1.d-10) then

        er(k)=1.d-10

       endif

      enddo

c  Lawley's approximation => a

      call phi(cov,a,v,e,er,w,n,nvec,io)

c  v=>s

      do j=1,nvec

       do i=1,n

        s(i,j)=v(i,j)

       enddo

      enddo

c  Objective function Sum(e-log(e)-1)

      z1=0.d0

      do k=kp1,n

       z1=z1+e(k)-dlog(e(k))-1.d0

      enddo

c  sqrt(er(e-1))*v

      do i=1,nvec

       do j=1,n

        ch=e(i)-1.d0

        if(ch<=0.d-10) then

         ch=1.d-10

        endif

        v(j,i)=s(j,i)*dsqrt(er(j)*ch)

       enddo

      enddo

      if(sw==1) then

c  Output eigenvalues

       call wrtfac(e,z1,nr,n,nvec,ni,io)

      endif

c-------

      deallocate(indx,w)

      return

      end

c

c===========================================

c

      subroutine wrtfac(e,z1,nr,n,nvec,ni,io)

c*************************************************

c  n - Dimension of sample covariance, nvec - Number of common factors

c  Eigenvalues of (er**-1/2)*cov*(er**-1/2) => e(n)

c=================================================

      real*8 e(n),ch,z1

      intent(in) e,z1,nr,n,nvec,ni,io

c*************************************************

c  Output eigenvectors

      write(io,'(/7x,''Maximum Likelihood Factor Analysis Eigenvalues'')

     *')

      write(io,'(/11x,''Eigenvalue'',5x,''%'',9x,''% Total''/)')

      t=0.

      sum=0.

      do k=1,n

       sum=sum+sngl(e(k))

      enddo

      do k=1,n

       s=100.*sngl(e(k))/sum

       t=t+s

       write(io,'(1x,i2,2x,f14.2,4x,2(f6.2,7x))') k,e(k),s,t

      enddo

      z2=dfloat(2*n+5)/6.d0+dfloat(2*nvec)/3.d0

      nd=((n-nvec)*(n-nvec)-(n+nvec))/2

      ch=z1*(dfloat(nr)-z2)

      write(io,'(/2x,''Chi Square Statistic = '',f12.1,2x,''Degrees of F

     *reedom =  '',i4)') sngl(ch),nd

      write(io,'(7x,''Number of Iterations = '',i4)') ni

c-------

      return

      end

c

c===========================================

c

      subroutine phi(cov,a,v,e,er,w,n,nvec,io)

c*************************************************

c  nvec - Number of common factors, cov - Sample covariance matrix

c  er(n) - Specific variances

c  e(n) - Eigenvalues (er**-1/2)*cov*(er**-1/2), v - eigenvectors

c=================================================

      implicit real*8(a-h,o-z)

      real*8 a(n,n),cov(n,n),v(n,n),e(n),er(n),w(n)

      intent(in) n,cov,er,io

      intent(out) a,e

c*************************************************

c  a=(er**-1/2)*cov*(er**-1/2)

      do i=1,n

       do j=i,n

        a(j,i)=cov(j,i)/dsqrt(er(i)*er(j))

        a(i,j)=a(j,i)

       enddo

      enddo

c  Eigenvalues & eigenvectors of a

      call eigen(a,n,n,e,w,v,jev)

c  a=(er**-1)-(er**-1/2)*(vv')*(er**-1/2)

      do i=1,n

c  Trap negative e's

       if(e(i)<=0.) then

        e(i)=1.d0

       endif

       do j=i,n

        z=0.d0

        do k=1,nvec

         z=z+v(i,k)*v(j,k)

        enddo

        a(j,i)=-z/dsqrt(er(i)*er(j))

        a(i,j)=a(j,i)

       enddo

       a(i,i)=a(i,i)+1.d0/er(i)

      enddo

c  a**2

      do i=1,n

       do j=i,n

        a(j,i)=a(j,i)*a(j,i)

        a(i,j)=a(j,i)

       enddo

      enddo

c-------

      return

      end

c

c===========================================

c

      subroutine fd(cov,v,e,er,n,nvec,f)

c********************************************

c  1st derivatives f(i)=del(S)/del(er(i))

c============================================

      implicit real*8(a-h,o-z)

      real*8 cov(n,n),v(n,n),e(n),er(n),f(n)

      intent(in) cov,v,e,er,n,nvec

      intent(out) f

c********************************************

      do i=1,n

      z=0.d0

       do j=1,nvec

        z=z+(e(j)-1.d0)*(v(i,j)*v(i,j))

       enddo

       f(i)=(z+1.d0-cov(i,i)/er(i))/er(i)

      enddo

c-------

      return

      end

c

c===========================================

c

      subroutine sd(v,e,er,n,nvec,f,a,d)

c********************************************

c  2nd derivative matrix of S => d

c============================================

      implicit real*8(a-h,o-z)

      real*8 v(n,n),e(n),er(n),f(n),a(n,n),d(n,n)

      integer r

      intent(in) a,e,er,f,v,n,nvec

      intent(out) d

c********************************************

      nc=nvec+1

      do i=1,n

       do j=i,n

        z=0.d0

        do k=1,nvec

         y=0.d0

         do r=nc,n

          y=y+2.d0*e(k)*(e(r)-1.d0)*v(i,r)*v(j,r)/(e(r)-e(k))

         enddo

         z=z+y*v(i,k)*v(j,k)

        enddo

        d(j,i)=a(j,i)+z

        d(i,j)=d(j,i)

       enddo

       d(i,i)=d(i,i)-2.d0*f(i)/er(i)

      enddo

c-------

      return

      end

c

c===========================================

c

      subroutine dludcmp(a,n,np,indx,d,ier,io)

c*************************************************

c  Computes the LU decomposition of the matrix a

c=================================================

      implicit real*8(a-h,o-z)

      parameter(nmax=500,tiny=1.0d-20)

      real*8 a(np,np),vv(nmax)

      real d

      integer indx(n),nmax

      intent(in) n,np,io

      intent(out) d,ier

c*************************************************

      d=1.

      ier=0

      do i=1,n

       aamax=0.d0

       do j=1,n

        if(dabs(a(i,j))>aamax) then

         aamax=dabs(a(i,j))

        endif

       enddo

       if(aamax==0.d0) then

        write(io,*) 'Singular matrix in ludcmp'

        ier=1

        return

       endif

       vv(i)=1.d0/aamax

      enddo

      do j=1,n

       do i=1,j-1

        sum=a(i,j)

        do k=1,i-1

         sum=sum-a(i,k)*a(k,j)

        enddo

        a(i,j)=sum

       enddo

       aamax=0.d0

       do i=j,n

        sum=a(i,j)

        do k=1,j-1

         sum=sum-a(i,k)*a(k,j)

        enddo

        a(i,j)=sum

        dum=vv(i)*dabs(sum)

        if(dum>=aamax) then

         imax=i

         aamax=dum

        endif

       enddo

       if(j/=imax) then

        do k=1,n

         dum=a(imax,k)

         a(imax,k)=a(j,k)

         a(j,k)=dum

        enddo

        d=-d

        vv(imax)=vv(j)

       endif

       indx(j)=imax

       if(dabs(a(j,j))<tiny) then

        a(j,j)=tiny

       endif

       if(j/=n) then

        dum=1.d0/a(j,j)

        do i=j+1,n

         a(i,j)=a(i,j)*dum

        enddo

       endif

      enddo

c--------

      return

      end

c

c===========================================

c

      subroutine dlubksb(a,n,np,indx,b)

c*************************************************

c  Solves the set of n linear equations Ax=b

c=================================================

      implicit real*8(a-h,o-z)

      integer indx(n),i,ii,j,ll

      real*8 a(np,np),b(n)

      intent(in) a,n,np

      intent(in out) b

c*************************************************

      ii=0

      do i=1,n

       ll=indx(i)

       sum=b(ll)

       b(ll)=b(i)

       if(ii/=0) then

        do j=ii,i-1

         sum=sum-a(i,j)*b(j)

        enddo

        elseif(sum/=0.) then

         ii=i

       endif

       b(i)=sum

      enddo

      do i=n,1,-1

       sum=b(i)

       do j=i+1,n

        sum=sum-a(i,j)*b(j)

       enddo

       b(i)=sum/a(i,i)

      enddo

c--------

      return

      end

c

c===========================================

c

      subroutine v2(m,v,cr,cns,r,ro,b1,b2,io)

c********************************************

c   Rotate candidate positions in two dimensions

c   ro - Rotation angle if uncontrained, no output if sw=0

c   m - Number of candidates used, v - Factor loadings for candidates

c  Output: cr - Rotated candidate positions

c============================================

      integer m,cns

      real*8 v(m-1,2)

      real cr(m,2)

      intent(in) m,v,cns,r,ro,b1,b2,io

      intent(out) cr

c********************************************

      if(cns/=2) then

       cosr=cos(ro)

       sinr=sin(ro)

        else

         cosr=cos(r)

         sinr=sin(r)

      endif

c  Candidate positions

      do k=1,m-1

c  Rotation

       s1=sngl(v(k,1))*cosr-sngl(v(k,2))*sinr

       s2=sngl(v(k,1))*sinr+sngl(v(k,2))*cosr

c  Scale by sigmas

       cr(k,1)=s1/b1

       cr(k,2)=s2/b2

      enddo

c  Candidate m

      cr(m,1)=0.

      cr(m,2)=0.

c-------

      return

      end

c

c===========================================

c

      subroutine param(b,cns,b1,b2,cosr,sinr,r,ro,xm,io)

c********************************************

c   Parameters from least squares fit of mean score

c   ro - rotation angle in radians, cosr - cosine(ro)

c   b1 & b2 - sigmas of x1 & x2, xm - mean ideal point

c   r - fixed rotation angle if cns=2, xm=0 if cns=1

c============================================

      real*8 b(5),del,s1,s2

      real xm(2),r

      integer cns

      intent(in) b,cns,r,io

      intent(out) b1,b2,cosr,sinr,ro,xm

c********************************************

      pi=3.141592

      s1=b(2)-b(1)

      s2=b(1)+b(2)

      if(dabs(s1)>1.d-14) then

       del=datan2(b(3),s1)/2.d0

       ro=sngl(del)

       del=dsin(2.d0*del)

       b1=sngl(s2-b(3)/del)/2.

       b2=sngl(s2+b(3)/del)/2.

       if(cns/=2) then

c  Rotation cosine & sin

        cosr=cos(sngl(del))

        sinr=sin(sngl(del))

       endif

        else

         if(cns/=2) then

          cosr=1./sqrt(2.)

          sinr=cosr

          ro=pi/2.

          b1=sngl(s2-b(3))/2.

          b2=sngl(s2+b(3))/2.

         endif

      endif

      if(cns==2) then

c  Rotation set to angle r

       angle=pi*r/180.

       ro=angle

       cosr=cos(angle)

       sinr=sin(angle)

      endif

c-------

c  Trap negative b's

      if(b1<=1.e-7) then

       write(io,*) 'Warning'

       write(io,*) 'b1 = ',b1,' < 0. b1 set to 1.e-4'

       b1=1.d-4

      endif

      if(b2<=1.e-7) then

       write(io,*) 'Warning'

       write(io,*) 'b2 = ',b2,' < 0. b2 set to 1.e-4'

       b2=1.d-4

      endif

c  Mean ideal points

      x1b=(b(4)*cosr-b(5)*sinr)/2.

      x2b=(-b(4)*sinr-b(5)*cosr)/2.

      b1=1./sqrt(b1)

      b2=1./sqrt(b2)

      if(cns==0) then

       xm(1)=x1b*b1

       xm(2)=x2b*b2

      endif

      if(cns==1) then

c  Mean=0

       xm(1)=0.

       xm(2)=0.

      endif

c-------

      return

      end

c

c===========================================

c

      subroutine regres(y,x,b,a,c,n,p,int,r2,sw,ier,io)

c*************************************************

c  Input: x - matrix of independent variables, y - dependent variables

c         n - No. of observations, p - No. of independent variables

c         int - Fit a mean if int = 1, sw=2 - output stats, sw=1 & 2 - header

c  Output: b - Regression coefficients, c - standard errors, r2

c=================================================

      implicit real*8(a-h,o-z)

      allocatable::s,w

      integer p,n,int,sw

      real r2

      real*8 y(n),x(n,p),b(p),a(p,p),c(p),w(:),s(:)

      intent(in) y,x,n,p,int,sw,io

      intent(out) b,c,r2,ier

c*************************************************

      allocate(s(0:p*p-1),w(0:p-1),stat=ibad)

      if(ibad/=0) then

       write(io,*) 'Unable to allocate work space for s in regres'

       return

      endif

      sy=0.d0

      do k=0,p-1

       w(k)=0.d0

      enddo

      if(int==1) then

       rn=dfloat(n)

c  Intercept

       do i=1,n

        sy=sy+y(i)

        do k=1,p

         w(k-1)=w(k-1)+x(i,k)

        enddo

       enddo

       sy=sy/rn

       do k=0,p-1

c  Mean of x(i,k)

        w(k)=w(k)/rn

       enddo

      endif

c-------

      do j=1,p

       do k=j,p

        a(j,k)=0.d0

       enddo

      enddo

c  Covariance matrix

      do i=1,n

       do j=1,p

        dij=x(i,j)-w(j-1)

        a(j,j)=a(j,j)+dij*dij

        do k=j+1,p

         dij=x(i,j)-w(j-1)

         dik=x(i,k)-w(j-1)

         a(j,k)=a(j,k)+dij*dik

         a(k,j)=a(j,k)

        enddo

       enddo

      enddo

      do k=1,p

       st=0.d0

       do i=1,n

        st=st+(x(i,k)-w(k-1))*(y(i)-sy)

       enddo

       c(k)=st

      enddo

      ssy=0.d0

      do i=1,n

       ssy=ssy+(y(i)-sy)*(y(i)-sy)

      enddo

c-------

c  a=>s

      do i=1,p

       do j=1,p

        st=a(j,i)

        s((i-1)*p+j-1)=st

       enddo

      enddo

      call dls(s,b,c,w,p,1,ier,io)

      if(ier==1) then

       write(io,*)'Covariance matrix in regres is not positive definite'

       return

      endif

      ca=0.d0

      if(int==1) then

c  Estimate intercept

       do k=1,p

        ca=ca+w(k-1)*b(k)

       enddo

      endif

      cr=sy-ca

      r=0.d0

      do i=1,n

       st=0.d0

       do j=1,p

        st=st+x(i,j)*b(j)

       enddo

       res=y(i)-st

       r=r+res*res

      enddo

c-------

      sdy=dsqrt(ssy/dfloat(n-1))

      r2=sngl(1.d0-r/ssy)

      se=dsqrt(r/dfloat(n-p))

      if(sw==1.or.sw==2) then

c  Output statistics

       write(io,'(2x,''Standard Deviation of Y ='',f15.3)') sdy

       write(io,'(2x,''R Square = '',f7.5,2x,''Standard Error = '',f15.4

     */)') r2,se

      endif

      if(int==1) then

       write(io,'(7x,''Intercept = '',g10.3/)') cr

      endif

c  Standard errors of OLS estimates

      do k=1,p

       c(k)=se*dsqrt(s((k-1)*p+k-1))

      enddo

      if(sw==2) then

       write(io,'(7x,''OLS Estimates'',6x,''Standard Error''/)')

       write(io,'(3x,i2,2x,g10.3,10x,g10.3)') (k,b(k),c(k),k=1,p)

       write(io,'(1x)')

      endif

c-------

      deallocate(s,w)

      return

      end

c

c===========================================

c

      subroutine strap(a,c,cov,d,dbar,e,er,f,s,v,w,m,mn,n,nu,ntry,vi,va,

     *cns,dim,bmn,br,bs1,bs2,bv,nstr,io)

c********************************************

c  Bootstap mean, rotation, variances using data in c

c  bmn - means & sigmas of mean, br - rotation, bs1 & bs2 - sigmas of x1 & x2

c  w - OLS estimates (b in regres & param)

c============================================

      allocatable::y

      real c(n,m),va(m),y(:,:),bmn(4),br(2),bs1(2),bs2(2),bv(2),xm(2)

      real*8 a(m,m),d(mn,mn),dbar(m),cov(mn,mn),e(m),er(m),f(mn),s(mn,mn

     *),v(mn,mn),w(mn*mn),s1,s2,s3,s4,rbt,z1

      integer cns,dim,sw,vi

      intent(in) c,dbar,va,m,mn,n,nu,nstr,ntry,vi,io

      intent(out) bmn,br,bs1,bs2,bv

c********************************************

      allocate(y(0:n-1,0:m-1),stat=ibad)

      if(ibad/=0) then

       write(io,*) 'Unable to allocate work space for y in strap'

       stop

      endif

      s1=0.d0

      s2=0.d0

      s3=0.d0

      s4=0.d0

      do k=1,2

       br(k)=0.

       bs1(k)=0.

       bs2(k)=0.

       bv(k)=0.

      enddo

c  Bootstrap c=>y

      wr=0

      kount=0

      do kb=1,nstr

       call boot(c,y,e,m,n,nu,io)

       call covm(y,cov,e,n,nu,mn,io)

c  Factor loadings =>s from maxle

       call maxle(cov,n,mn,dim,a,d,e,er,f,s,w,sge,z1,ntry,0,ier,io)

       if(ier==0) then

        kount=kount+1

c  Set variables for OLS fit & solve for parameters

        do k=1,mn

         cov(k,1)=s(k,1)*s(k,1)

         cov(k,2)=s(k,2)*s(k,2)

         cov(k,3)=s(k,1)*s(k,2)

         cov(k,4)=s(k,1)

         cov(k,5)=s(k,2)

        enddo

        ir=5

c  Check for valence weight

        if(vi>0) then

         ir=6

         do k=1,mn

          cov(k,6)=va(k)

         enddo

        endif

c  OLS fit dbar=cov*w

         int=0

         if(kb==nstr) then

c  Write r*2

          sw=1

           else

            sw=0

         endif

         call regres(dbar,cov,w,d,er,mn,ir,int,r2,sw,ier,io)

c  Solve for parameters from w

         call param(w,cns,b1,b2,cosr,sinr,r,ro,xm,io)

c  Sum parameter estimates

         s1=s1+dble(xm(1))

         s2=s2+dble(xm(2))

         s3=s3+dble(xm(1))*dble(xm(1))

         s4=s4+dble(xm(2))*dble(xm(2))

         br(1)=br(1)+ro

         br(2)=br(2)+ro*ro

         bs1(1)=bs1(1)+b1

         bs1(2)=bs1(2)+b1*b1

         bs2(1)=bs2(1)+b2

         bs2(2)=bs2(2)+b2*b2

         if(vi>0) then

          bv(1)=bv(1)+w(6)

          bv(2)=bv(2)+w(6)*w(6)

         endif

       endif

      enddo

c  End of bootstrap loop

      if(nstr<2) then

       bmn(1)=sngl(s1)

       bmn(2)=sngl(s2)

       return

      endif

      if(kount==0) then

       write(io,'(7x,/''ERROR - each bootstrap has an error''/)')

       stop

      endif

      rbt=dfloat(kount)

c  Mean & standard deviations

      bmn(1)=sngl(s1/rbt)

      bmn(2)=sngl(s2/rbt)

      bmn(3)=sngl(dsqrt(s3/rbt-s1/rbt*s1/rbt))

      bmn(4)=sngl(dsqrt(s4/rbt-s2/rbt*s2/rbt))

      br(1)=br(1)/rbt

      br(2)=br(2)/rbt-br(1)*br(1)

      br(2)=sqrt(br(2))

      bs1(1)=bs1(1)/rbt

      bs1(2)=bs1(2)/rbt-bs1(1)*bs1(1)

      bs1(2)=sqrt(bs1(2))

      bs2(1)=bs2(1)/rbt

      bs2(2)=bs2(2)/rbt-bs2(1)*bs2(1)

      bs2(2)=sqrt(bs2(2))

      if(vi>0) then

       bv(1)=bv(1)/rbt

       bv(2)=bv(2)/rbt-bv(1)*bv(1)

       bv(2)=sqrt(bv(2))

      endif

      write(io,'(7x,''Number of Boostraps = '',i3)') kount

c-------

      deallocate(y)

      return

      end

c

c===========================================

c

      subroutine boot(c,y,e,m,n,nu,io)

c********************************************

c  Bootstapped x signal in y

c============================================

      allocatable::w

      real c(n,m),y(n,m),w(:)

      real*8 e(m)

      intent(in) c,m,n,nu

      intent(out) e,y

c********************************************

      allocate(w(0:nu-1),stat=ibad)

      if(ibad/=0) then

       write(io,*) 'Unable to allocate work space for w in boot'

       stop

      endif

      rn=float(nu)

c  Uniforms

      call random_number(w)

      do k=1,m

       e(k)=0.d0

      enddo

      do k=1,nu

       nsam=max(int(w(k-1)*rn),1)

       do j=1,m

        y(k,j)=c(nsam,j)

        e(j)=e(j)+dble(y(k,j))

       enddo

      enddo

c  Mean score differences

      do k=1,m-1

       e(k)=e(k)/dfloat(nu)

      enddo

c------

      deallocate(w)

      return

      end

c

c===========================================

c

      subroutine v1(m,v,cr,z,d,e,b,er,f,vi,va,wv,cns,bs,ier,io)

c********************************************

c   Estimate candidate positions in one dimension

c    m - Number of candidates used, v - Factor loadings

c    vi= 0 - No candidate characteristic values used

c      = 1 - Candidate characteristic weight estimated

c    va - Candidate characteristic values derived from survey data

c    cns = all parameters estimated in OLS fit

c         = 1 mean ideal point & the reference candidate = 0

c  Output: cr - candidate positions, wv - valence weight, bs - sig(x)

c=================================================

      real*8 v(m-1),e(m),b(m),d(m-1,m-1),z(m-1,6),er(m),f(*)

      real cr(m,1),va(m),r2,wv

      integer cns,vi

      intent(in) m,e,v,vi,va,cns,io

      intent(out) bs,cr,wv,ier

c*************************************************

      mn=m-1

      ir=2

c  Set variables for OLS fit

      do k=1,mn

       z(k,1)=v(k)*v(k)

       z(k,2)=v(k)

      enddo

c  Check for valience weight

      if(vi>0) then

       ir=ir+1

       do i=1,mn

        z(i,ir)=va(i)

       enddo

      endif

c-------

c  OLS fit - no intercept

      call regres(e,z,b,f,er,mn,ir,1,r2,0,ier,io)

c  Solve for parameters

      bs=1./b(1)

      x1b=0.

      if(cns==0) then

       x1b=-0.5/b(2)

      endif

      if(bs<=1.e-7) then

       write(io,*) 'Warning'

       write(io,*) 'bs = ',bs,' < 0. Absolute value is used'

      endif

      bs=sqrt(abs(bs))

      x1b=x1b*bs

c  Estimate of candidate positions

      do k=1,m

       cr(k,1)=v(k)/bs-x1b

      enddo

      if(vi>0) then

       wv=b(ir)

        else

         wv=0.

      endif

      cr(m,1)=-x1b

c-------

      return

      end

c

c===========================================

c

      subroutine eigen(cov,np,n,d,e,z,jev)

c********************************************

c  Finds eigenvalues and eigenvectors of a symmetric real matrix COV

c   whose order is n. np is the row dimension of cov set in the

c   program dimension (or parameter) statement

c  Output:

c    d contains the eigenvalues in descending order. If jev = j, then

c     eigenvalues are correct but unordered for indices 1,2,..,jev-1.

c    jth column of z is the j-th orthonormal eigenvector.

c    jev = zero  for normal return,

c        = j     if the j-th eigenvalue has not been determined after

c                 30 iterations.

c============================================

      real*8 cov(np,n),d(n),e(n),z(np,n)

      intent(in) cov,np,n

      intent(out) d,z,jev

c********************************************

      if(n<1) then

       write(6,*) 'Matrix order n < 2 in subroutine eigen'

       return

        elseif(n>np) then

         write(6,*) 'Matrix order n > limit',np,'in subroutine eigen'

         return

      endif

c-------

c  Call tridiagonalizing subroutine

      call tridag(cov,d,e,z,n,np)

c  Call main eigenvalue program

      call evs(d,e,z,n,np,jev)

c--------

      return

      end

c

c============================================

c

      subroutine tridag(a,d,e,z,n,np)

c********************************************

c  This subroutine reduces a real symmetric matrix to a symmetric

c  tridiagonal matrix using orthogonal similarity transformations.

c  Input:

c     n is the order of the matrix;

c     a contains the real*8 symmetric input matrix;

c  Output:

c     d contains the diagonal elements of the tridiagonal matrix;

c     e contains the subdiagonal elements of the tridiagonal

c       matrix in its last n-1 positions. e(1) is set to zero;

c     z contains the orthogonal transformation matrix

c       produced in the reduction;

c============================================

      integer i,j,k,l,n,ii,jp1

      real*8 a(np,n),d(n),e(n),z(np,n),f,g,h,hh,scale

c********************************************

      do i=1, n

       d(i)=0.d0

       do j = 1, i

        z(i,j)=a(i,j)

       enddo

      enddo

      do 30 ii = 2, n

         i = n + 2 - ii

         l = i - 1

         h = 0.d0

         scale = 0.d0

         if (l < 2) go to 17

c     :::::::::: scale row ::::::::::

        do 12 k = 1, l

  12    scale = scale + dabs(z(i,k))

c

        if (scale /= 0.d0) go to 14

  17    e(i) = z(i,l)

         go to 290

c

  14    do 15 k = 1, l

            z(i,k) = z(i,k) / scale

            h = h + z(i,k) * z(i,k)

  15    continue

c

         f = z(i,l)

         g = -dsign(dsqrt(h),f)

         e(i) = scale * g

         h = h - f * g

         z(i,l) = f - g

         f = 0.d0

c

         do 24 j = 1, l

            z(j,i) = z(i,j) / h

            g = 0.d0

c     :::::::::: form element of a*u ::::::::::

            do 180 k = 1, j

  180       g = g + z(j,k) * z(i,k)

c

            jp1 = j + 1

            if (l < jp1) go to 220

c

            do 20 k = jp1, l

  20       g = g + z(k,j) * z(i,k)

c     :::::::::: form element of p ::::::::::

  220       e(j) = g / h

            f = f + e(j) * z(i,j)

  24  continue

c

         hh = f / (h + h)

c     :::::::::: form reduced a ::::::::::

         do 260 j = 1, l

            f = z(i,j)

            g = e(j) - hh * f

            e(j) = g

c

            do 260 k = 1, j

               z(j,k) = z(j,k) - f * e(k) - g * z(i,k)

  260    continue

c

  290    d(i) = h

  30  continue

c     :::::::::: accumulation of transformation matrices ::::::::::

      do 500 i = 1, n

         l = i - 1

         if (d(i) == 0.d0) go to 380

c

         do 360 j = 1, l

            g = 0.d0

c

            do 340 k = 1, l

  340       g = g + z(i,k) * z(k,j)

c

            do 360 k = 1, l

               z(k,j) = z(k,j) - g * z(k,i)

  360    continue

c

  380    d(i) = z(i,i)

         z(i,i) = 1.d0

         if (l < 1) go to 500

c

         do 400 j = 1, l

            z(i,j) = 0.d0

            z(j,i) = 0.d0

  400    continue

c

  500 continue

c-------

      return

      end

c

c============================================

c

      subroutine evs(d,e,z,n,np,jev)

c********************************************

c  Input:

c   n is the order of the matrix;

c   d contains the diagonal elements of the input matrix;

c   e contains the subdiagonal elements of the input matrix

c     in its last n-1 positions. e(1) is arbitrary;

c   z contains the transformation matrix produced in the

c     reduction by TRIDAG if called. If the eigenvectors

c     of the tridiagonal matrix are desired, z must contain

c     the identity matrix.

c  Output:

c   d contains the eigenvalues in descending order. If an

c     error exit is made, the eigenvalues are correct but

c     unordered for indices 1,2,...,jev-1;

c   e has been destroyed;

c   z contains orthonormal eigenvectors of the symmetric

c     tridiagonal (or full) matrix. If an error exit is made,

c     z contains the eigenvectors associated with the stored ev's;

c   jev = zero for a normal return

c       = j  if the j-th eigenvalue has not been determined

c            after 30 iterations.

c============================================

      integer i,j,k,l,m,n,ii,mml,jev

      real*8 d(n),e(n),z(np,n),b,c,f,g,p,r,s,machep

c---------------------------------

c     machep is a machine dependent parameter specifying

c     the relative precision of floating point arithmetic.

c     machep = 16.0d0**(-13) for long form arithmetic on s360

      data machep/16.d-40/

c********************************************

      jev = 0

      do 10 i = 2, n

  10  e(i-1) = e(i)

c

      e(n) = 0.d0

c

      do 24 l = 1, n

       j = 0

c     :::::::::: look for small sub-diagonal element ::::::::::

  1   do 11 m = l, n

        if (m == n) go to 12

        if (dabs(e(m))<=machep * (dabs(d(m))+dabs(d(m+1)))) goto 12

c

  11    continue

c

  12   p = d(l)

        if (m == l) goto 24

c     set error -- no convergence to an eigenvalue after 30 iterations

      if (j == 30) then

          jev=1

          return

      endif

         j = j + 1

c     :::::::::: form shift ::::::::::

         g = (d(l+1) - p) / (2.d0 * e(l))

         r = dsqrt(g*g+1.d0)

         g = d(m) - p + e(l) / (g + dsign(r,g))

         s = 1.d0

         c = 1.d0

         p = 0.d0

         mml = m - l

c     :::::::::: for i=m-1 step -1 until l do -- ::::::::::

         do 20 ii = 1, mml

          i = m - ii

          f = s * e(i)

          b = c * e(i)

           if (dabs(f) < dabs(g)) goto 15

          c = g / f

          r = dsqrt(c*c+1.d0)

          e(i+1) = f * r

          s = 1.d0 / r

          c = c * s

          goto 16

  15      s = f / g

          r = dsqrt(s*s+1.d0)

          e(i+1) = g * r

          c = 1.d0 / r

          s = s * c

  16      g = d(i+1) - p

          r = (d(i) - g) * s + 2.d0 * c * b

          p = s * r

          d(i+1) = g + p

          g = c * r - b

c     :::::::::: form vector ::::::::::

          do 18 k = 1, n

           f = z(k,i+1)

           z(k,i+1) = s * z(k,i) + c * f

           z(k,i) = c * z(k,i) - s * f

  18       continue

c

  20  continue

c

         d(l) = d(l) - p

         e(l) = g

         e(m) = 0.d0

         go to 1

  24  continue

c     :::::::::: order eigenvalues and eigenvectors ::::::::::

      do 30  ii = 2, n

       i = ii - 1

       k = i

       p = d(i)

c

       do 26  j = ii, n

       if (d(j) > p) then

        k = j

        p = d(j)

       endif

  26   continue

c

         if (k /= i) then

          d(k) = d(i)

          d(i) = p

c

          do 28 j = 1, n

           p = z(j,i)

           z(j,i) = z(j,k)

           z(j,k) = p

  28     continue

        endif

c

  30  continue

c-------

      return

      end

c

c============================================

c

      subroutine dls(a,x,b,d,n,inv,ier,io)

c********************************************

c  Solves the n linear equations Ax=b

c  Input: a - positive definite matrix

c  Output: x - solution, a - inverse of a matrix if inv = 1

c============================================

      real*8 a(n,n),b(n),d(n),x(n),sum

      intent(in) b,n,inv,io

      intent(out) x,ier

      intent(in out) a

c********************************************

      ier=0

      do i=1,n

       do j=i,n

        sum=a(i,j)

        do k=i-1,1,-1

         sum=sum-a(i,k)*a(j,k)

        enddo

        if(i==j) then

         if(sum<=0.d0) then

          ier=1

c         write(io,*) 'a is not positive definite in DLS'

          return

         endif

         d(i)=dsqrt(sum)

          else

           a(j,i)=sum/d(i)

        endif

       enddo

      enddo

c-------

      do i=1,n

       sum=b(i)

       do k=i-1,1,-1

        sum=sum-a(i,k)*x(k)

       enddo

       x(i)=sum/d(i)

      enddo

      do i=n,1,-1

       sum=x(i)

       do k=i+1,n

        sum=sum-a(k,i)*x(k)

       enddo

       x(i)=sum/d(i)

      enddo

      if(inv/=1) then

       return

      endif

c-------

c  Lower diagonal inverse of L in a

      do i=1,n

       a(i,i)=1.d0/d(i)

       do j=i+1,n

        sum=0.d0

        do k=i,j-1

         sum=sum-a(j,k)*a(k,i)

        enddo

        a(j,i)=sum/d(j)

       enddo

      enddo

c-------

c  Inverse of LL'

      do i=1,n

       do j=i+1,n

        sum=0.d0

        do k=j,n

         sum=sum+a(k,i)*a(k,j)

        enddo

        a(i,j)=sum

       enddo

c--------

c Diagonal

       sum=0.d0

       do k=i,n

        sum=sum+a(k,i)*a(k,i)

       enddo

        a(i,i)=sum

      enddo

      do i=1,n

       do j=i+1,n

        a(j,i)=a(i,j)

       enddo

      enddo

c-------

      return

      end

c

c===========================================

c

      subroutine ideal2(n,nu,m,x,id,typ,uv,c,dbar,va,wv,a,e,w,y,flg,chk,

     *vd,dav,pow,ref,sc,ni,pr,io)

c********************************************

c  Fit ideal points in a two dimensional space. Version 3-16-2004

c  x - Input: Transformed scores selected in dat (c in call). Best is 0

c  typ - type code for respondents to be plotted

c  flg - flag for typ code, n - Number of respondents, m - Candidates

c  vd - respondent id's

c  uv - flag for candidate whose votes are to be estimated

c  Output: Ideal points. x(i,3) = kt if one or more candidates have a 0 score,

c  x(i,3) = 11 if the line segment fit is used, x(i,3) = 21 if the least

c  squares fit is used

c  ni - number used, a - Sum of voters nearest to each candidate

c  c - Matrix of candidate positions

c  chk - magnitude of difference between max & min x score

c  xmn =>y - mean of estimated ideal points

c============================================

      real x(n,m),c(m,2),va(m),y(7),wv,sc

      real*8 a(m),dbar(m),e(m)

      integer id(nu),typ(n),uv(m),w(m),flg(7),pow,pr

      character*5 vd(nu)

      character*2 ref

      character dav

      intent(in) n,nu,m,typ,c,dbar,va,wv,flg,chk,vd,dav,pow,ref,sc,uv,pr

     *,io

      intent(out) a,id,ni,y

      intent(in out) x

c*************************************************

      mn=m-1

      do k=1,5

       y(k)=0.

      enddo

c  X'X of candidate points

      do k=1,mn

       smx=c(k,1)-c(m,1)

       smn=c(k,2)-c(m,2)

       y(1)=y(1)+smx*smx

       y(2)=y(2)+smn*smn

       y(3)=y(3)+smx*smn

      enddo

c  Determinant*2

      rm=(y(1)*y(2)-y(3)*y(3))*2.

      ni=0

c-------

c  Initialize seed

      call random_seed

      xm1=0.

      xm2=0.

c  Find minimum & maximum of the ith row transformed scores (smn,smx)

      do i=1,nu

       nxi=typ(i)

       if(nxi==flg(1).or.nxi==flg(2).or.nxi==flg(3).or.nxi==flg(4).or.nx

     *i==flg(5).or.nxi==flg(6).or.nxi==flg(7).or.flg(1)==-1) then

        ni=ni+1

        id(ni)=i

c  Add back last column

        do j=1,mn

         e(j)=x(i,j)+x(i,m)

        enddo

        e(m)=x(i,m)

        smn=e(1)

        smx=smn

        kt=0

        do j=1,m

         st=e(j)

         if(st<smn) then

          smn=st

         endif

         if(st>smx) then

          smx=st

         endif

c-------

c  Find j's for zero transformed scores

         if(st==0.) then

c  Count zero scores

          kt=kt+1

          w(kt)=j

         endif

        enddo

c  Set ideal point at the candidate with a 0. If 1<kt then set the ideal

c   point at the mean of the candidate positions.

        call random_number(y(6:7))

        if(kt==1) then

         x(id(ni),1)=c(w(kt),1)+sc*(y(6)-.5)

         x(id(ni),2)=c(w(kt),2)+sc*(y(7)-.5)

         x(id(ni),3)=float(kt)

          elseif(kt>1) then

           x1=0.

           x2=0.

           do k=1,kt

            x1=x1+c(w(k),1)

            x2=x2+c(w(k),2)

           enddo

           x(id(ni),1)=x1/float(kt)+sc*(y(6)-.5)

           x(id(ni),2)=x2/float(kt)+sc*(y(7)-.5)

           x(id(ni),3)=float(min(kt,10))

        endif

c--------

c  Convex combination on line

        kmn=0

        kmx=0

        if(kt==0) then

         do j=1,m

          as=e(j)

          if(abs(as-smn)<1.e-7) then

c  j at min

           kmn=kmn+1

c  j has the highest score but not max value

           w(1)=j

            elseif(abs(as-smx)<1.e-7) then

c  j at max

             kmx=kmx+1

c  j has the lowest score at the min value of transformed grades

             w(2)=j

          endif

         enddo

c-------

         if(kmn==1 .and. kmx==1 .and. smx>=chk) then

          if(pow==2) then

           a1=sqrt(e(w(1)))

           a2=sqrt(e(w(2)))

            else

             a1=e(w(1))

             a2=e(w(2))

          endif

c  Ratio of inverse grades

          st=-(a2/a1)

          dc=(c(w(1),1)-c(w(2),1))**2+(c(w(1),2)-c(w(2),2))**2

          ar=st/sqrt(dc)

c  Prorate on cord between thetas

          call random_number(y(6:7))

          x(id(ni),1)=ar*c(w(2),1)+(1.-ar)*c(w(1),1)+sc*(y(6)-.5)

          x(id(ni),2)=ar*c(w(2),2)+(1.-ar)*c(w(1),2)+sc*(y(7)-.5)

          x(id(ni),3)=11.

           else

c  Least squares fit using the adjusted scores

            y(4)=0.

            y(5)=0.

            do k=1,mn

             e(k)=x(i,k)-dbar(k)

             y(4)=y(4)+(c(k,1)-c(m,1))*e(k)

             y(5)=y(5)+(c(k,2)-c(m,2))*e(k)

            enddo

c  Ideal points

            x(id(ni),1)=(y(3)*y(5)-y(2)*y(4))/rm+c(m,1)

            x(id(ni),2)=(y(3)*y(4)-y(1)*y(5))/rm+c(m,2)

            x(id(ni),3)=21.

         endif

        endif

       endif

      enddo

c-------

      write(io,'(2x,70(''*''))')

      if(flg(1)>-1) then

       write(io,'(/7x,i5,'' respondent ideal points were selected using

     *the flags:''/)') ni

       write(io,'(7('' Plot('',i2'')= '',i2))') (k,flg(k),k=1,7)

       if(ni==0) then

        stop

       endif

        else

        write(io,'(/7x,''All respondent ideal points are estimated''/)')

      endif

      do k=1,m

       a(k)=0.

      enddo

c  Find  minimum distance from the kth voter to the candidates

      do k=1,ni

       smn=1.e7

       knc=0

       do nc=1,m

        if(uv(nc)==1) then

c  Distance from x(k) to c(nc)

         distc=(c(nc,1)-x(id(k),1))**2+(c(nc,2)-x(id(k),2))**2+wv*va(nc)

         if(distc<0.) then

c  Dist=0

          dist=0.

         endif

         if(distc>=0.) then

          dist=sqrt(distc)

         endif

         if(dist<smn) then

          smn=dist

          knc=nc

         endif

        endif

       enddo

c  Sum for candidate with smallest distance from i's ideal point

       a(knc)=a(knc)+1.

       ncode=nint(x(id(k),3))

       if(ref(1:1)=='y') then

c  Reflect around axis 1

        x1=-1.*x(id(k),1)

         else

          x1=x(id(k),1)

       endif

       if(ref(2:2)=='y') then

c  Reflect around axis 2

        x2=-1.*x(id(k),2)

         else

          x2=x(id(k),2)

       endif

c  Sum x1 & x2

       xm1=xm1+x1

       xm2=xm2+x2

c-------

       if(pr==1) then

        if(dav=='y') then

c  Output ideal points to file 7

          write(7,'(a5,1x,2(f10.2,1x),1x,i2)') vd(id(k)),x1,x2,ncode

           else

            write(7,'(i5,1x,2(f10.2,1x),1x,i2)') id(k),x1,x2,ncode

        endif

       endif

      enddo

c  End of min distance loop for new method

      xm1=xm1/float(ni)

      xm2=xm2/float(ni)

      write(io,'(7x,''Mean of the Estimated Ideal Points = ('',f12.2,'',

     *'',f12.2,'')''/)') xm1,xm2

      write(io,'(2x,70(''=''))')

      y(1)=xm1

      y(2)=xm2

c-------

      return

      end

c

c===========================================

c

      subroutine ideal1(n,nu,m,x,c,typ,w,sc,flg,ni,io)

c********************************************

c  Fit ideal points in a one dimensional space for candidate positions c

c  x - Input: Scores  Output: Ideal points (Note: array c is cr in call)

c  nu - Number of respondents use, m - Number of candidates

c============================================

      real x(n,m),c(m),smx,smn,st,d,sc,y

      integer typ(n),w(m),flg(5)

      intent(in) n,m,c,sc,io

      intent(in out) x

c*************************************************

c Initialize seed

      call random_seed

      ni=0

c  Find minimum & maximum of the ith row scores

      do i=1,nu

       nti=typ(i)

c  Fit ideal point for i if flags match or flg(1)=-1

       if(nti==flg(1).or.nti==flg(2).or.nti==flg(3).or.nti==flg(4).or.nt

     *i==flg(5)) then

        ni=ni+1

c  Add back last column

        do j=1,m-1

         x(i,j)=x(i,j)+x(i,m)

        enddo

c-------

        smn=x(i,1)

        smx=smn

        kt=0

        do j=1,m

         st=x(i,j)

         if(st<smn) then

          smn=st

         endif

         if(st>smx) then

          smx=st

         endif

c--------

c  Find j's for zero scores

         if(st==0.) then

          kt=kt+1

          w(kt)=j

         endif

        enddo

c  Set ideal point at the candidate with a 0. If 1<kt<4 then set the

c   ideal point at the mean of the candidate positions.

        call random_number(y)

        if(kt==1) then

         x(i,1)=c(w(kt))+sc*y

         x(i,2)=float(kt)

          elseif(kt>1) then

           x1=0.

           do k=1,kt

            x1=x1+c(w(k))

           enddo

           x(i,1)=x1/float(kt)+sc*y

           x(i,2)=float(min(kt,10))

            elseif(kt==0) then

             do j=1,m

              st=x(i,j)

              if(abs(st-smn)<1.e-7) then

               w(1)=j

                elseif(abs(st-smx)<1.e-7) then

                 w(2)=j

              endif

             enddo

c-------

c  Ideal point from quadratic model

             d=c(w(1))-c(w(2))

             st=(x(i,w(1))-x(i,w(2)))/(2.*d*d)

             x(i,1)=(.5-st)*c(w(1))+(.5+st)*c(w(2))+sc*y

             x(i,2)=11.

        endif

       endif

      enddo

c-------

      return

      end

c

c===========================================

c

      function charead(par,string,delim,ia,ier,io)

c*******************************************

c  Find a parameter name in the character*700 string. Then find the

c  character*50 string charead in string enclosed by the symbols in delim.

c  There may be blanks before delim(1:1) and up to two blanks after it.

c  ier=1 - parameter not found in string, ier=2 - no left delimiter,

c  ier=3 - no right delimiter

c===========================================

      character*50 charead

      character*700 string

      character*20 par

      character*2 delim

      intent(in) par,string,delim,io

      intent(out) ia,ier

c*******************************************

      ier=0

c  Find par

      nc=len_trim(par)

      ib=index(string,par(1:nc))

      if(ib==0) then

       ier=1

       return

      endif

      charead='          '

      ia=index(string(ib:),delim(1:1))

      if(ia==0) then

       ier=2

       return

      endif

c  Position of first character in par

      ia=ia+ib

      ic=index(string(ia:),delim(2:2))-1

      if(ic<=0) then

       ier=3

       return

      endif

      charead=string(ia:ia+ic)

c----------

      return

      end

c

c===========================================

c

      function rnumread(par,string,delim,ia,ier,io)

c*******************************************

c  Find a parameter name in the character*700 string. Then find the real

c  number in the string enclosed by the symbols in delim.

c  There may be blanks before delim(1:1) and up to two blanks after it.

c  ier=1 - parameter not found in string, ier=2 - no left delimiter,

c  ier=3 - no right delimiter, ier=4 - no period, ier=5 - number > 999,999

c===========================================

      real rnumread

      real*8 mag,frac

      character*700 string

      character*20 par

      character*10 ac,buf

      character*2 delim

      intent(in) par,delim,io

      intent(out) ia,ier

c*******************************************

      ier=0

c  Find par

      nc=len_trim(par)

      ib=index(string,par(1:nc))

      if(ib==0) then

       ier=1

       return

      endif

      ia=index(string(ib:),delim(1:1))

      if(ia==0) then

       ier=2

       return

      endif

c  Position of first character in par

      ia=ia+ib

c  Trap blanks after delim(1:1)

      ic=index(string(ia:ia),' ')

      if(ic>0) then

       ia=ia+1

      endif

      ic=index(string(ia:ia),' ')

      if(ic>0) then

       ia=ia+1

      endif

c  Check for negative

      if(string(ia:ia)=='-') then

       rneg=-1.

       ia=ia+1

        else

         rneg=1.

      endif

      it=index(string(ia:),'.')

      if(it==0) then

       ier=4

       return

      endif

c  Size of number

      ic=it-1

      if(ic<8) then

       call number(string(ia:),ic,numag,io)

        else

         write(io,*) 'The number ',string(ia:),' > 9,999,999'

         ier=5

         return

      endif

      ib=ia+it

      buf=string(ib-1:)

c  Check on blank after period

      if(buf(2:2)==' ') then

       buf(2:2)='0'

      endif

      ic=index(buf(2:),delim(2:2))-1

      if(ic<=0) then

       ier=3

       return

      endif

      if(ic<7) then

       ac=buf(2:)

       call number(ac,ic,numfrac,io)

        else

         write(io,*) 'The number ',string(ib:),' > 9,999,999'

         ier=5

         return

      endif

      mag=dfloat(numag)

      frac=dfloat(numfrac)/dfloat(10**ic)

      rnumread=sngl(mag+frac)*rneg

c----------

      return

      end

c

c=============================================

c

      subroutine number(string,sn,num,io)

c*******************************************

c  Integer from character integer in string at pointer ip

c  sn - field for integer num

c===========================================

      integer sn,num

      character*7 string

      intent(in) string,sn,io

      intent(out) num

c********************************************

      select case(sn)

       case(1)

        num=ichar(string(1:1))-48

       case(2)

        n1=ichar(string(1:1))-48

        n2=ichar(string(2:2))-48

        num=n1*10+n2

       case(3)

        n1=ichar(string(1:1))-48

        n2=ichar(string(2:2))-48

        n3=ichar(string(3:3))-48

        num=n1*100+n2*10+n3

       case(4)

        n1=ichar(string(1:1))-48

        n2=ichar(string(2:2))-48

        n3=ichar(string(3:3))-48

        n4=ichar(string(4:4))-48

        num=n1*1000+n2*100+n3*10+n4

       case(5)

        n1=ichar(string(1:1))-48

        n2=ichar(string(2:2))-48

        n3=ichar(string(3:3))-48

        n4=ichar(string(4:4))-48

        n5=ichar(string(5:5))-48

        num=n1*10000+n2*1000+n3*100+n4*10+n5

       case(6)

        n1=ichar(string(1:1))-48

        n2=ichar(string(2:2))-48

        n3=ichar(string(3:3))-48

        n4=ichar(string(4:4))-48

        n5=ichar(string(5:5))-48

        n6=ichar(string(6:6))-48

        num=n1*100000+n2*10000+n3*1000+n4*100+n5*10+n6

       case(7)

        n1=ichar(string(1:1))-48

        n2=ichar(string(2:2))-48

        n3=ichar(string(3:3))-48

        n4=ichar(string(4:4))-48

        n5=ichar(string(5:5))-48

        n6=ichar(string(6:6))-48

        n7=ichar(string(7:7))-48

        num=n1*1000000+n2*100000+n3*10000+n4*1000+n5*100+n6*10+n7

      endselect

c-------

      return

      end

c

c===========================================

c

      function numread(par,string,delim,ia,ier,io)

c*******************************************

c  Find a parameter name in the character*700 string. Then find the

c  integer numread in string enclosed by the symbols in delim.

c  There may be blanks before delim(1:1) and up to two blanks after it.

c  ier=1 - parameter not found in string, ier=2 - no left delimiter,

c  ier=3 - no right delimiter, ier=4 - number > 9999999

c===========================================

      integer numread

      character*700 string

      character*20 par

      character*2 delim

      intent(in) par,string,delim,io

      intent(out) ia,ier

c*******************************************

      ier=0

c  Find par

      nc=len_trim(par)

      ib=index(string,par(1:nc))

      if(ib==0) then

       ier=1

       return

      endif

      ia=index(string(ib:),delim(1:1))

      if(ia==0) then

       ier=2

       return

      endif

c  Position of first character in par

      ia=ia+ib

c  Trap blanks after delim(1:1)

      ic=index(string(ia:ia),' ')

      if(ic>0) then

       ia=ia+1

      endif

      ic=index(string(ia:ia),' ')

      if(ic>0) then

       ia=ia+1

      endif

c  Check for negative

      if(string(ia:ia)=='-') then

       neg=-1

       ia=ia+1

        else

         neg=1

      endif

c  Size of number

      ic=index(string(ia:),delim(2:2))-1

      if(ic<=0) then

       ier=3

       return

      endif

      if(ic<8) then

       call number(string(ia:),ic,numread,io)

        numread=numread*neg

        else

         write(io,*) 'The number ',string(ia:),' > 999999'

         ier=4

         return

      endif

c----------

      return

      end

c

c=============================================

c

      subroutine flagread(string,delim,flg,nc,ier,io)

c*******************************************

c  Find integer flags in string

c===========================================

      integer flg(nc)

      character*50 string

      character*10 par

      character*2 delim

      intent(in) string,delim,io

      intent(out) flg,ier

c*******************************************

      ier=0

      do k=1,nc

       if(k<10) then

        par=achar(k+48)

         else

          mn=k/10

          par(1:1)=achar(mn+48)

          par(2:2)=achar(k-mn*10+48)

       endif

       na=len_trim(par)

       ib=index(string,par(1:na))

       if(ib==0) then

        ier=1

        return

       endif

       ia=index(string(ib:),delim(1:1))

       if(ia==0) then

        ier=2

        return

       endif

c  Position of first character in par

       ia=ia+ib

c  Trap blanks after delim(1:1)

       ic=index(string(ia:ia),' ')

       if(ic>0) then

        ia=ia+1

       endif

       ic=index(string(ia:ia),' ')

       if(ic>0) then

        ia=ia+1

       endif

c  Check for negative

       if(string(ia:ia)=='-') then

        neg=-1

        ia=ia+1

         else

          neg=1

       endif

c  Size of number

       ic=index(string(ia:),delim(2:2))-1

       if(ic<=0) then

        ier=3

        return

       endif

       if(ic<7) then

        call number(string(ia:),ic,numread,io)

        nf=numread*neg

         else

          write(io,*) 'The number ',string(ia:),' > 999999'

          ier=4

          return

       endif

       flg(k)=nf

      enddo

c-------

      return

      end

c

c=============================================

c

      subroutine lineread(buf,names,mu,io)

c********************************************

c  Read names in buf

c============================================

      character*700 buf

      character*30 names(mu)

      intent(in) buf,mu,io

      intent(out) names

c********************************************

c  Check for character in first space

      nc=ichar(buf(1:1))

      if(nc/=32) then

       write(io,'(17x,''ERROR''/)')

       write(io,'(7x,''Leave at least one space before the first name in

     * the list of '',i2,'' candidate names'',/,6x,''above the first dat

     *a line'')') mu

       write(io,'(/(a80))') buf

       return

      endif

      npt=1

      do nkt=1,mu

c  Search for first non blank

       kb=npt

       dowhile(ichar(buf(kb:kb))==32)

        kb=kb+1

       enddo

c  Search for first blank after nkt name

       nc=kb

       dowhile(ichar(buf(nc:nc))/=32)

        nc=nc+1

       enddo

       names(nkt)=buf(kb:nc-1)

       npt=nc

      enddo

c-------

      return

      end