!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!    PROGRAM TEST98
!
!    AUTHORS: Valeria Ruggiero (Department of Mathematics - University of Ferrara)
!             Luca Zanni (Department of Mathematics - University of Modena)
!             Carla Durazzi (Department of Mathematics - University of Ferrara)
!
!    This program generates a sparse quadratic programming (QP) problem
!    with linear constraints of the form
!
!                      minimize    f(x)=1/2 xt G x + qt x
!                      subject to  A x >=d1
!                                  C x =d2
! 
!    where x is the n-dimensional unknown vector;
!    m is the number of all the contraints; G is an n*n matrix;
!    A is an (m-neq)*n matrix; A is an neq*n matrix of full row-rank and 
!         _   _        _  _
!        |  A  |      | d1 |   
!    B = |     |   d= |    |
!        |_ C _|      |_d2_|
!        
!    From now on, K() will be used to denote the spectral condition number 
!    of a matrix and Z is the matrix whose columns are an orthogonal base 
!    for the null space of the matrix of the active constraints.
!    In input, the user has to introduce the data described below (see  
!    "Description of the input variables" and file inputtest that is 
!    an example of data input).
!     
!    The FUNCTION RANF for the generation of random numbers is processor 
!    dependent: its substitution with a suitable one is left to the user.
!
!    This program uses double precision floating point arithmetic: in the description
!    the word "real" stands for double precision real.
!
!--------------- DESCRIPTION OF THE INPUT VARIABLES ----------------------------------
!
! iseme            = seed for the generation of random numbers               (integer)
! n                = n. of variables (order of the matrix G)                 (integer)
! m                = n. of all the constraints                               (integer)
! nac              = n. of active constraints + n. of equality constraints   (<= m, n; integer)
! neq              = n. of equality constraints                              (<=n; integer)
! condg            = log_10 (K(G)))                                          (integer)  
! condb            = log_10 (K(B))                                           (integer)
! condzgz          = log_10 (K(Zt G Z))                                      (integer)
! condba           = log_10 (K(matrix of the active constraints))            (integer)
! glmin,ngrango    = minimum non-zero eigenvalue and rank of matrix G        (real,integer)
! zgzlmin,nzgzrango= minimum non-zero eigenvalue and rank of matrix Zt*G*Z   (real,integer)
! blmin            = minimum non-zero singular value of matrix B             (real)
! balmin           = minimum non-zero singular value of the matrix of the active 
!                    constraints                                             (real)
! indg, indb       = two different distributions for the singular values of the matrices G 
!                    and B are allowed: if indg (or indb) is equal to zero a logarithmic 
!                    distribution is used, if indg (or indb) is equal to one the uniform one,
!                    is employed, otherwise they are equally spaced between the smallest 
!                    and the largest singular value (see subroutine SING; integer,integer).
! ndeg             = level of degeneracy of the constraints                  (integer). 
!                    In fact the multipliers are computed by the following rule            
!                               lambda(i)=10**(-ranf()*ndeg)       
! sparg, sparb     = sparsity of the matrices G and B (percentage)           (real,real)
! nome             = name for the file containing the generated matrix G     (character)
! nomevin          = name for the file containing the generated constraints  (character)
! indequi          = if indequi is non-zero, indequi is the bandwidth of 
!                    the matrices G and B; if it is equal to zero, no control is given
!                    when generating the matrices.                           (integer)
!
!------------------------------------------------------------------------------
!---------- OUTPUT: DESCRIPTION -----------------------------------------------
!
! As output, the program gives two unformatted files which will be described below. 
! First we point out the way the sparse matrices G and B are stored: each row contains
! only the non-zero elements of G (or B), then a vector kg (or kb) is used to keep the 
! number of non-zero elements of each row and a matrix ig (or ib) is arranged so that 
! ig(i,*) contains the column index of the elements in the i-th row of G.
! EXAMPLE:
!        _        _             _  _        _       _      _       _
!       |  3  0  2 |           |  2 |      | 1  3    |    | 3  2    | 
!    G= |  0  4  1 |  ==>  kg= |  2 |   ig=| 2  3    |  g=| 4  1    | 
!       |_ 2  1  4_|           |_ 3_|      |_1  2  3_|    |_2  1  4_| 
!
! Two unformatted files are produced as output of this program. 
! The first one (=nome) contains :
! - n,ic (=n. of all the non-zero elements of the upper triangular part of G) and the
!   rank of G                          (integer)
! - the matrix G stored in a set of consecutive records which contain the non-zero
!   elements of the upper triangular part of G, in a row wise order. Each record has the 
!   informations of four non-zero elements of G in the following format
!       (row index i, column index j, element (i,j) of the matrix)  
! - the vector q                       (real)
! - the minimum value of the objective function f(x) (real)
! - the sparsity of G                  (integer)
! The second file (=nomevin) contains:
! - m,neq,nac-neq,ic (=n. of non-zero elements of the matrix of the constraints)
! - all the non-zero elements of the matrix B in the same order as those of G.  
! - the solution x of the problem      (real)
! - the vectors b and d                (real)
! - the sparsity of B                  (integer)
!    
!    
!******* REFERENCES *********************************************************************
!
!    1) Ruggiero V., Zanni L.; "On a class of Iterative Methods for Large-
!       Scale  Convex Quadratic Programs",Rend. del Circ. Matem. di Palermo, Serie II,
!       Suppl. 58 (1999), pp. 213-227.
!    2) Stewart G.W.; "The efficient Generation of Random Orthogonal Matrices
!       with an Application to Condition Estimators",SIAM J. Numer. Anal., 17 
!       n.3 (1980), pp. 403-409.
!
!*****************************************************************************************
!
!                                 08/04/98
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      program TEST98
      parameter(nmax=7001,igmax=400,mmax=7001,ibmax=400,
     +           ndim=10001,!max(nmax,mmax)
     +           kmax=2,lmax=20000,ntot=ndim*igmax)
	implicit double precision (a-h,o-z)

      real*8  b(mmax,ibmax),g(nmax,igmax),vcos(lmax),
     +      x(0:nmax),lambda(0:mmax),q(nmax),d(mmax)
      integer ib(mmax,ibmax),ig(nmax,igmax),kb(mmax),kg(nmax),
     +         irv(lmax),jcv(lmax)
      character*15 nome,nomevin


      common/work/fcoeff(ntot),irf(ntot),icf(ntot)
	common /seme/iseme

!--------- EXECUTABLE STATEMENTS ------------------------------------------

 
      tol=1d-15

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------- DATA INPUT --------------------------------------------------
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
!       write(*,*)'Give the seed for the generation of random numbers'
      read(*,*) iseme
!       write(*,*)'Introduce the dimension n of the problem and '
!     1           ,'the total number of the constraints'
        read(*,*) n,m
!40   write(*,*)'Introduce the number of equalities (neq <=n) and the '
!     1           ,'number of the active constraints (nac <= m,n)'
40    read(*,*) nac,neq
      if ((neq.gt.n).or.(nac.gt.min(m,n))) then
	write(*,*)'neq <=n and nac<=m,n'
	goto 40
      endif


!       write(*,*)'log_10 (K(G))= ,log_10 (K(B))= '
      read(*,*)condg,condb
!       write(*,*)'log_10 (K(ZT*G*Z))'
!     1,',log_10 (K(matrix of the active constraints))'
      read(*,*)condzgz,condba
!       write(*,*)'Introduce the minimum non-zero eigenvalue of the '
!     1,'matrices G, B, Zt G Z'
!       write(*,*)'and of the  matrix of the active constraints'
      read(*,*)glmin,blmin,zgzlmin,balmin
!       write(*,*)'If you want a  logarithmic distribution for elements '
!     1,'in G, press 0,'
!       write(*,*)'otherwise an uniform distribution is used' 
       read(*,*)indg
!       write(*,*)'If you want a  logarithmic distribution for elements in B, press 0,'
!       write(*,*)'otherwise an uniform distribution is used' 
       read(*,*)indb
!       write(*,*)'Introduce the rank of G and Zt*G*Z'
       read(*,*)ngrango,nzgzrango
!       write(*,*)'The level of degeneracy is ....'
       read(*,*)ndeg
!       write(*,*)'Introduce the sparsity of the matrices G and B (percentage)'
       read(*,*)sparg,sparb
!       write(*,*)'The name of the file containing the generated matrix G is ....'
       read(*,'(a15)')nome                
!       write(*,*)'The name of the file containing the constraints is ....'
       read(*,'(a15)')nomevin             
!       write(*,*)'Introduce the number of non-zero elements for each row of the randomly' 
!       write(*,*)'generated matrices G and B; if it is equal to zero, no control is given '
!       write(*,*)'when generating these matrices'
       read(*,*) indequi 
       if(indequi.eq.0)indequi=n

       nulld1=n-nac-nzgzrango             ! zero diagonal elements of the matrix D1
       nulld2=n-ngrango-nulld1            ! zero diagonal elements of the matrix D2
       nnullb=M*n*(1.-sparb/100.)         ! non-zero diagonal elements of the matrix B
       nnullg=n*n*(1.-sparg/100.)         ! non-zero diagonal elements of the matrix G

       write(*,*)'This program generates the following problem'
       write(*,*)
       write(*,*)'minimize    1/2 xt G x + qt x '
       write(*,*)'subject to  C x >=d'                      
       write(*,*)'            A x =b'
       write(*,*)   
       write(*,*)   
       write(*,*)'where ',n,' is the order of G and '
       write(*,*)m,' is the number of all the constraints'
       write(*,*)'n. of active constraints= ',nac-neq
       write(*,*)'n. of equalities= ',neq
       write(*,*)'K(G)= 10**',condg,' rank(G)= ',ngrango
       write(*,*)' min. eigenvalue of G= ',glmin
       write(*,*)
       write(*,*)'     _   _ '
       write(*,*)'    |  A  |'
       write(*,*)'B = |     |'
       write(*,*)'    |_ C _|'
       write(*,*)
       write(*,*)'K(B)= 10**',condb,'min. eigenvalue of B= ',blmin
       write(*,*)'K(Zt G Z) = 10**',condzgz,' rank(Zt G Z)= ',
     1         nzgzrango,' min. eigenvalue of Zt G Z= ',zgzlmin
       write(*,*)
       if(indg.eq.1)then
	 write(*,*)'Uniform distribution for the eigenvalues of G'
       else if(indg.eq.0)then
	 write(*,*)'Logarithmic distribution for the eigenvalues of G'
       else
	 write(*,*)'The eigenvalues of G are equally spaced '
     1,'between the smallest and the largest one'
       endif
!
       if(indb.eq.1)then
	 write(*,*)'Uniform distribution for the singular values of B'
       else if(indb.eq.0)then
	 write(*,*)'Logarithmic distribution for the singular '
     1    ,'values of B'
       else
	 write(*,*)'The singular values of G are equally spaced '
     1    ,'between the smallest and the largest one'
       endif
       write(*,*)
       write(*,*)'Level of degeneracy= ',ndeg
       write(*,*)
       write(*,*)'Sparsity of G =',sparg,' Sparsity of B =',sparb
       write(*,*)       
       write(*,*)'The matrix G will be stored in: ',nome
       write(*,*)'The matrix B will be stored in: ',nomevin
       write(*,*) 
!
!
!---------- DATA INPUT: END -----------------------------------------
!
!
!
!----- Generating G --------------------------------------------------
!
	 ieqg=nnullg/n  +1

       do i=1,n
	 do j=1,igmax
	 g(i,j)=0.0
	 ig(i,j) =0
	 enddo
       enddo
       icontg=ngrango
       icontdg=ngrango
       do i=1,n
	 kg(i)=1
	 ig(i,1)=i
       enddo

       call sing(g,nac,nac-nulld2,condg,glmin,indg)
       if(n-nac.ne.0)then
          call sing(g(nac+1,1),n-nac,nzgzrango,condzgz,zgzlmin,indg)
       endif
       k=0
1007   if(icontdg.lt.n)then
	 do ir=1,n/2
	   if(abs(g(ir,1)).lt.tol)then
	     jc=ranf()*(n-ir)+1+ir
	     goto 345
	   endif
	 enddo
	 do jc=n/2+1,n
	   if(abs(g(jc,1)).lt.tol)then
	     ir=ranf()*(jc-1)+1
	     goto 345
	   endif
	 enddo
       endif
       if(icontg.gt.nnullg)goto 2001
234    ii=0
2344   ir=ranf()*(n-n/10)+1
       jc=ranf()*(min(n-ir,indequi))+1+ir
       ii=ii+1
	  if ((kg(ir).gt.ieqg.or.kg(jc).gt.ieqg)
     1          .and.ii.lt.1000)goto 2344
	  if(ii.ge.1000)then
            write(*,*)'sparg is too big'
            stop
          endif
        if(ir.eq.jc)goto 2344

345    k=k+1 
       if (k.gt.lmax) then
	  write(*,*)' error: lmax too small'  
	       stop
       endif   
       irv(k)=ir
       jcv(k)=jc
       vc=ranf()*2.-1.
       vcos(k)=vc
       call tleft(g,ig,kg,n,vc,icontg,icontdg,tol,nmax,igmax,ir,jc)
       call trigth(g,ig,kg,n,vc,icontg,icontdg,tol,nmax,igmax,ir,jc)

	 do j=2,kg(ir)
	  ii=ig(ir,j)
	  do l=1,kg(ii)
	     if(ig(ii,l).eq.ir)then
	       g(ii,l)=g(ir,j)
	       goto 333
	     endif
	  enddo
333       continue
	 enddo

	 do j=2,kg(jc)
	  ii=ig(jc,j)
	  do l=1,kg(ii)
	     if(ig(ii,l).eq.jc)then
	       g(ii,l)=g(jc,j)
	       goto 3334
	     endif
	  enddo
3334    continue
	 enddo

	 goto 1007

2001   write(*,*)'Non-zero elements of G=',icontg
	 do i=1,n
	  x(i)=kg(i)
	 enddo
	 write(*,*)
	 write(*,*)'maximum number of non-zero elements per row '
       write(*,*)'of G (kg) =',x(idamax(n,x(1),1))
	 write(*,*)
	 write(*,*)'Non-zero diagonal elements of G=',icontdg
	 write(*,*)'Obtained sparsity=',(1.-icontg/real(n*n))*100
       sparg=(1.-icontg/real(n*n))*100
	 write(*,*)
!
!----- Generating G: END --------------------------------------------
!
!
!----- Generating B -------------------------------------------------
!
      
       kb(1:m)=0
       do i=1,m
	 do j=1,ibmax
	   b(i,j)=0.0
	   ib(i,j)=0
	 enddo
       enddo
       mn=min(m,n)
       icontb=mn  
       icontdb=mn 
       do i=1,mn   
	 kb(i)=1
	 ib(i,1)=i
       enddo

       call sing(b,nac,nac,condba,balmin,indb)
	 if (m.gt.nac) then
       call sing(b(nac+1,1),m-nac,mn-nac,condb,blmin,indb)
	 endif

       ndtot=m*n*(1.-sparb/100.)
	   iequi=ndtot/m +1       

       do i=1,k
       call tleft(b,ib,kb,m,vcos(i),icontb,icontdb,tol,
     +            mmax,ibmax,irv(i),jcv(i))
       enddo

       do jc=n+1,m
121     ii=0
1211      ir=ranf()*(jc-nac-1) + nac+1
	  ii=ii+1
	  if (kb(ir).gt.iequi.and. ii.lt.1000) goto 1211
          if(ii.eq.1000)then
             write(*,*)'sparb id too big'
             stop
         endif
       ucos=ranf()*2.-1.
       call trigth(b,ib,kb,n,ucos,icontb,icontdb,tol,mmax,ibmax,ir,jc)
       enddo
       if(icontb.ge.ndtot)goto 100

999    continue 
	  ii=0
1212      ir=ranf()*(m-m/10)+1
	   ii=ii+1
	   if (kb(ir).gt.iequi.and. ii.lt.1000) goto 1212
           if(ii.eq.1000)then
                write(*,*)'sparb is too big'
           endif
	 if (ir.lt.nac) then
	  jc=ranf()*(min(indequi, nac-ir))+1+ir
	 elseif (ir.eq.nac)then
	     goto 999
	 else
	  jc=ranf()*(min(m-ir,indequi))+1+ir
	 endif

       ucos=ranf()*2.-1.

       call trigth(b,ib,kb,n,ucos,icontb,icontdb,tol,mmax,ibmax,ir,jc)
       if(icontb.ge.ndtot)goto 100
       goto 999

100    write(*,*)'n. of elements of  B =',icontb
       write(*,*)'n. trasf V=',k
       write(*,*)'non-zero diagonal elements of B=',icontdb
       write(*,*)'Obtained sparsity=',(1.-icontb/real(m*n))*100
       sparb=(1.-icontb/real(m*n))*100
       do i=1,m
	 d(i)=kb(i)
       enddo
       write(*,*)
       write(*,*)'maximum number of non-zero elements per row '
       write(*,*)'of B=',d(idamax(m,d,1))
       write(*,*)
!
!----- Generating B: END --------------------------------------------------
!
!----- Generating the solution x and the non-zero multiplier lambda: ------
!
       do i=1,n
	 x(i)=ranf()*2.-1.       
       enddo


       do i=1,nac
	 lambda(i)=10**(-ranf()*ndeg)
       enddo
!
!----- Generating the solution x and the non-zero multiplier lambda: END ---
!
!----------------- Generating d and b --------------------------------------
!
      
       d(1:m)=0.0
       do i=1,m
	 do j=1,kb(i)
	  d(i)=d(i)+b(i,j)*x(ib(i,j))
	 enddo
       enddo
       do i=nac+1,m
	 d(i)=d(i)-10*ranf()
       enddo
!
!----- Generating d and b: END -----------------------------------------
!
!----- Generating the vector q ----------------------------------
!
	
	 q(1:n)=0.0
	 do i=1,n
	 do j=1,kg(i)
	   q(i)=q(i)-g(i,j)*x(ig(i,j))
	 enddo
	 enddo
	 fmin=0.
	 do i=1,n
	 fmin=fmin-q(i)*x(i)
	 enddo
	 fmin=fmin*0.5
	 do i=1,nac
	 do j=1,kb(i)
	 q(ib(i,j))=q(ib(i,j))+lambda(i)*b(i,j)
	 enddo
	 enddo

!
!----- Generating the vector q: END --------------------------------------
!
!----- Compute the minimum value of the function -------------------------        
	 do i=1,n
	 fmin=fmin+q(i)*x(i)
	 enddo
!-----------------------------------------------------------------------        
!        
!---------- Writing the output files -----------------------------------
	 ic=0
	 do i=1,n
	 do j=1,kg(i)
	  if(ig(i,j).ge.i)then
	    ic=ic+1
	    irf(ic)= i
	    icf(ic)=ig(i,j)
	    fcoeff(ic)=g(i,j)
	  endif
	 enddo
	 enddo
	 open(unit=2,file=nome,form='unformatted')
	 write(2)n,ic,ngrango
	 nn=ic/4
	 do i=1,nn
	  write(2) (irf(j),icf(j),fcoeff(j),j=(i-1)*4+1,i*4)
	 enddo
       irest=ic-nn*4
       if (irest .ne. 0) then
	 write(2) (irf(j),icf(j),fcoeff(j),j=nn*4+1,ic)
       endif
       write(2)(q(i),i=1,n)
       write(2)fmin 
       write(2) sparg
	 close(2)
	 ic=0
	 do i=neq+1,m
	 do j=1,kb(i)
	 ic=ic+1
	 irf(ic)= i-neq
	 icf(ic)=ib(i,j)
	 fcoeff(ic)=b(i,j)
	 enddo
	 enddo
	 do i=1,neq   
	 do j=1,kb(i)
	 ic=ic+1
	 irf(ic)= m-neq+i
	 icf(ic)=ib(i,j)
	 fcoeff(ic)=b(i,j)
	 enddo
	 enddo
	 open(unit=7,file=nomevin,form='unformatted')
	 write(7)m,neq,nac-neq,ic
	 nn=ic/4
	 do i=1,nn
	 write(7) (irf(j),icf(j),fcoeff(j),j=(i-1)*4+1,i*4)
	 enddo
       irest=ic-nn*4
       if (irest .ne. 0) then
	 write(7) (irf(j),icf(j),fcoeff(j),j=nn*4+1,ic)
       endif
	 write(7)(x(i),i=1,n)
	 write(7)(d(i),i=neq+1,m),(d(i),i=1,neq)
	 write(7)sparb
	 close(7)
	 stop
	 end
!
!================= END PROGRAM TEST98 =======================================
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                  SUBROUTINES
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!====================== Subroutine sing =====================================
! If the matrix is square, it computes its eigenvalues, otherwise its singular 
! values
!
       subroutine sing(a,n,nrango,conda,almin,ind)
       implicit double precision (a-h,o-z)
       real*8  a(n)
       integer nrango
	   common /seme/iseme
       do i=nrango+1,n
	  a(i)=0.0
       enddo
       a(1)=10**conda
       a(nrango)=1.

!--------  Logarithmic distribution of the singular values ------------------

       if (ind.eq.0) then
	 do i=2,nrango-1
	   a(i)=10**(real(nrango-i)/real(nrango-1)*conda)
	 enddo

!--------  Uniform distribution of the singular values ------------------

       elseif (ind.eq.1)then
	 do i=2,nrango-1
	    t=ranf()
	   a(i)=t*(a(1)-a(nrango))+1
	 enddo
	else

!--------  Equally spaced singular values (between the smallest and the largest one) ------------------

	   ra=(a(1)-a(nrango))/(nrango-1)
	 do i=2,nrango-1
	   a(i)=a(i-1)-ra
	 enddo
       endif
	 do i=1,nrango
	   a(i)=a(i)*almin
	 enddo
       return
       end

!============= End subroutine sing =====================================
!
!============ Subroutine tleft =========================================
! It makes a Givens' transformation on matrix A: 
!     A* Gt(ir,jc)
! using random numbers in [-1,1] as non-zero elements of G
!
       subroutine tleft(a,ia,kr,n,c,icont,icontd,tol,lda,imax,ir,jc)
       parameter(nmax=7001,igmax=400,mmax=7001,
     +           ibmax=400,ndim=10001,kmax=2,
     +           lmax=20000,ntot=ndim*igmax)
       implicit double precision (a-h,o-z)
       real*8  a(lda,imax)
       integer ia(lda,imax),kr(imax)
       common/work/tempi(ndim),tempj(ndim),dum(ndim,igmax-2),
     +       indi(ndim),indj(ndim),idum(ndim,2*igmax-2)
	   common /seme/iseme
!
!     n = number of rows of A
!
       c2=c**2
       s2=1.-c2
       s=sqrt(s2)

      
	 tempi(1:n)=0.
	 tempj(1:n)=0.
	 indi(1:n)=0
	 indj(1:n)=0
     
       do k=1,n
       do i=1,kr(k)
	 if(ia(k,i).eq.ir)then
	    tempi(k)=a(k,i)
	    indi(k)=i
	 endif
	 if(ia(k,i).eq.jc)then
	    tempj(k)=a(k,i)
	    indj(k)=i
	 endif
       enddo
       enddo
       do k=1,n
       if(k.ne.ir)then
	 coeff=c*tempi(k)-s*tempj(k)
	      if(dabs(coeff).gt.tol)then
	       if(indi(k).ne.0)then
	       icol=indi(k)
	       else
	       kr(k)=kr(k)+1
	       if (kr(k).ge.igmax) then
	       write(*,*)'igmax is too small'
	       stop
	       endif
	       icol=kr(k)
	       icont=icont+1
	       ia(k,icol)=ir
	       endif
	      a(k,icol)=coeff
	      else
	       if(indi(k).ne.0)then
	       icont=icont-1
	       j=indi(k)
		  do i=j,kr(k)-1
		   a(k,i)=a(k,i+1)
		   ia(k,i)=ia(k,i+1)
		  enddo
	       a(k,kr(k))=0.0
	       ia(k,kr(k))=0
	    if(indj(k).ne.0.and.indj(k).gt.indi(k))indj(k)=indj(k)-1
	       indi(k)=0
	       kr(k)=kr(k)-1
	       endif
	      endif
	 endif
	 if(k.ne.jc)then
	 coeff=s*tempi(k)+c*tempj(k)
	      if(abs(coeff).gt.tol)then
	       if(indj(k).ne.0)then
	       icol=indj(k)
	       else
	       kr(k )=kr(k )+1
	       if (kr(k ).ge.igmax) then
	       write(*,*)'igmax is too small'
	       stop
	       endif
	       icol=kr(k )
	       icont=icont+1
	       ia(k ,icol)=jc
	       endif
	      a(k ,icol)=coeff
	      else
	       if(indj(k).ne.0)then
	       icont=icont-1
	       j=indj(k)
	      do i=j,kr(k )-1
		   a(k ,i)=a(k ,i+1)
		   ia(k ,i)=ia(k ,i+1)
		  enddo
	       a(k ,kr(k ))=0.0
	       ia(k ,kr(k ))=0
	   if(indi(k).ne.0.and.indi(k).gt.indj(k))indi(k)=indi(k)-1
	       indj(k)=0
	       kr(k )=kr(k )-1
	       endif
	      endif
	 endif
	 enddo !k  
	 if (ir.le.n) then
	   coeff=c*tempi(ir)-s*tempj(ir)
	   if(abs(coeff).gt.tol)then
	   if (abs(a(ir,1)).le.tol) then
	   icont=icont+1
	   icontd=icontd+1
	   endif
	   a(ir,1)=coeff
	   else
	   if (abs(a(ir,1)).gt.tol) then
	   icont=icont-1
	   icontd=icontd-1
	   endif
	   a(ir,1)=0.0   
	   endif
	  endif
	  if (jc.le.n) then
	   coeff=s*tempi(jc)+c*tempj(jc)
	   if(abs(coeff).gt.tol)then
	   if (abs(a(jc,1)).le.tol) then
	   icont=icont+1
	   icontd=icontd+1
	   endif
	   a(jc,1)=coeff
	   else
	   if (abs(a(jc,1)).gt.tol) then
	   icont=icont-1
	   icontd=icontd-1
	   endif
	   a(jc,1)=0.0   
	   endif
	  endif
	   return
	   end
!============ End subroutine tleft =========================================
!
!============ Subroutine tright =========================================
! It makes a Givens' transformation on matrix A: 
!     G(ir,jc) * A
! using random numbers in [-1,1] as non-zero elements of G
!

       subroutine trigth(a,ia,kr,n,c,icont,icontd,tol,lda,imax,ir,jc)
       parameter(nmax=7001,igmax=400,mmax=7001,
     +           ibmax=400,ndim=10001,kmax=2,
     +           lmax=20000,ntot=ndim*igmax)
       implicit double precision (a-h,o-z)
       real*8 a(lda,imax)
       integer ia(lda,imax),kr(imax)
       common/work/tempi(ndim),tempj(ndim),dum(ndim,igmax-2),
     +       indi(ndim),indj(ndim),idum(ndim,2*igmax-2)
	   common /seme/iseme
!
!    n=number of columns of A
!
       c2=c**2
       s2=1.-c2
       s=sqrt(s2)

     
	 tempi(1:n)=0.
	 tempj(1:n)=0.
	 indi(1:n)=0
	 indj(1:n)=0
     
       do i=1,kr(ir)
       tempi(ia(ir,i))=a(ir,i)
       indi(ia(ir,i))=i
       enddo
       do i=1,kr(jc)
       tempj(ia(jc,i))=a(jc,i)
       indj(ia(jc,i))=i
       enddo
       do k=1,n
	 if(k.ne.ir)then
	 coeff=c*tempi(k)-s*tempj(k)
	      if(abs(coeff).gt.tol)then
	       if(indi(k).ne.0)then
	       icol=indi(k)
	       else
	       kr(ir)=kr(ir)+1
	       if (kr(ir).ge.igmax) then
	       write(*,*)'igmax is too small'
	       stop
	       endif
	       icol=kr(ir)
	       icont=icont+1
	       ia(ir,icol)=k
	       endif
	      a(ir,icol)=coeff
	      else
	       if(indi(k).ne.0)then
	       icont=icont-1
	       j=indi(k)
		  do i=j,kr(ir)-1
		   a(ir,i)=a(ir,i+1)
		   ia(ir,i)=ia(ir,i+1)
		  enddo
	       a(ir,kr(ir))=0.0
	       ia(ir,kr(ir))=0
	       do i=1,n
	       if(indi(i).gt.j)indi(i)=indi(i)-1
	       enddo
	       indi(k)=0
	       kr(ir)=kr(ir)-1
	       endif
	      endif
	 endif
	 if(k.ne.jc)then
	 coeff=s*tempi(k)+c*tempj(k)
	      if(abs(coeff).gt.tol)then
	       if(indj(k).ne.0)then
	       icol=indj(k)
	       else
	       kr(jc)=kr(jc)+1
	       if (kr(jc).ge.igmax) then
	       write(*,*)'igmax is too small'
	       stop
	       endif
	       icol=kr(jc)
	       icont=icont+1
	       ia(jc,icol)=k
	       endif
	      a(jc,icol)=coeff
	      else
	       if(indj(k).ne.0)then
	       icont=icont-1
	       j=indj(k)
		  do i=j,kr(jc)-1
		   a(jc,i)=a(jc,i+1)
		   ia(jc,i)=ia(jc,i+1)
		  enddo
	       a(jc,kr(jc))=0.0
	       ia(jc,kr(ir))=0
	       do i=1,n
	       if(indj(i).gt.j)indj(i)=indj(i)-1
	       enddo
	       indj(k)=0
	       kr(jc)=kr(jc)-1
	       endif
	      endif
	 endif
	 enddo !k  

	 if (ir.le.n) then
	   coeff=c*tempi(ir)-s*tempj(ir)
	   if(abs(coeff).gt.tol)then
	   if (abs(a(ir,1)).le.tol) then
	   icont=icont+1
	   icontd=icontd+1
	   endif
	   a(ir,1)=coeff
	   else
	   if (abs(a(ir,1)).gt.tol) then
	   icont=icont-1
	   icontd=icontd-1
	   endif
	   a(ir,1)=0.0   
	   endif
	  endif
	  if (jc.le.n) then
	   coeff=s*tempi(jc)+c*tempj(jc)
	   if(abs(coeff).gt.tol)then
	   if (abs(a(jc,1)).le.tol) then
	   icont=icont+1
	   icontd=icontd+1
	   endif
	   a(jc,1)=coeff
	   else
	   if (abs(a(jc,1)).gt.tol) then
	   icont=icont-1
	   icontd=icontd-1
	   endif
	   a(jc,1)=0.0   
	   endif
	  endif
	   return
	   end
!============ Subroutine tright =====================================================
!
!========= FUNCTION RANF ============================================================
! This function generates random numbers: it is processor dependent, because function
! RAN is processor dependent.

	   double precision function ranf()
	   common /seme/iseme
  10       ranf=ran(iseme)
	   if (ranf.eq.0) goto 10
	   return 
	   end
!=====================================================================================