subroutine cg_cheb(hoyre,alf,bet,flag,N,m1,m2)
  !
  !
  !      The routine solves the generalized biharmonic equation using
  !      Shen's special Chebyshev basis. The capacitance matrix, P, is
  !      solved using the Bi-CGSTAB method.
  !
  !
  use Interface, ONLY: cg_cheb_hele,amvprod,cheb_eig,cheb_sepfakt,  &
       cg_cheb_pfirst,cg_cheb_pre_pfirst,compDiag,  &
       leg_to_shen,matr_to_vekt,vekt_to_matr  
  !
  use Bi_Parameter
  !
  implicit none
  !
  integer, intent(in)          :: N,m1,m2,flag
  double precision, intent(in) :: alf
  double precision, intent(in) :: bet
  double precision, intent(inout) :: hoyre(0:N-1,0:N-1)
  !
  integer i,ii,q,j,odde,l
  integer flagget
  !
  integer, save :: oldn=0
  !
  double precision rk
  !
  !    for timing purposes only
  !
  integer count,ocount,count_rate,count_max

  real tt(2)
  !
  !
  !      All variables that are to be saved for subsequent 
  !      right hand sides.
  !
  !
  double precision, dimension(:),allocatable, save :: rr1
  double precision, dimension(:),allocatable, save :: rr2
  double precision, dimension(:),allocatable, save :: rr3
  double precision, dimension(:),allocatable, save :: ss2
  double precision, dimension(:),allocatable, save :: ss3
  double precision, dimension(:),allocatable, save :: D1
  double precision, dimension(:),allocatable, save :: p1
  double precision, dimension(:),allocatable, save :: q1
  double precision, dimension(:),allocatable, save :: r1
  double precision, dimension(:),allocatable, save :: s1
  double precision, dimension(:),allocatable, save :: D2
  double precision, dimension(:),allocatable, save :: p2
  double precision, dimension(:),allocatable, save :: q2
  double precision, dimension(:),allocatable, save :: r2
  double precision, dimension(:),allocatable, save :: s2
  double precision, dimension(:),allocatable, save :: dd
  double precision, dimension(:),allocatable, save :: e1
  double precision, dimension(:),allocatable, save :: e2
  double precision, dimension(:),allocatable, save :: e3
  double precision, dimension(:),allocatable, save :: e4
  double precision, dimension(:),allocatable, save :: wsave
  double precision, dimension(:),allocatable, save :: Sum1e1
  double precision, dimension(:),allocatable, save :: Sum2e1
  double precision, dimension(:),allocatable, save :: Sum3e1
  double precision, dimension(:),allocatable, save :: Sum4e1
  double precision, dimension(:),allocatable, save :: Sum1e2
  double precision, dimension(:),allocatable, save :: Sum2e2
  double precision, dimension(:),allocatable, save :: Sum3e2
  double precision, dimension(:),allocatable, save :: Sum4e2
  double precision, dimension(:),allocatable, save :: Sum1e3
  double precision, dimension(:),allocatable, save :: Sum2e3
  double precision, dimension(:),allocatable, save :: Sum3e3
  double precision, dimension(:),allocatable, save :: Sum4e3
  double precision, dimension(:),allocatable, save :: Sum1e4
  double precision, dimension(:),allocatable, save :: Sum2e4
  double precision, dimension(:),allocatable, save :: Sum3e4
  double precision, dimension(:),allocatable, save :: Sum4e4
  double precision, dimension(:),allocatable, save :: Diag1
  double precision, dimension(:),allocatable, save :: Diag2
  double precision, dimension(:),allocatable, save :: Dinv1
  double precision, dimension(:),allocatable, save :: Dinv2
  !
  double precision, dimension(:,:),allocatable, save :: Id1
  double precision, dimension(:,:),allocatable, save :: Id1inv
  double precision, dimension(:,:),allocatable, save :: Id2
  double precision, dimension(:,:),allocatable, save :: Id2inv
  double precision, dimension(:,:),allocatable, save :: Ja1
  double precision, dimension(:,:),allocatable, save :: Fa1
  double precision, dimension(:,:),allocatable, save :: Ja2
  double precision, dimension(:,:),allocatable, save :: Fa2
  double precision, dimension(:,:),allocatable, save :: Sum1
  double precision, dimension(:,:),allocatable, save :: Sum2
  double precision, dimension(:,:),allocatable, save :: Stor1Up
  double precision, dimension(:,:),allocatable, save :: Stor2Up
  double precision, dimension(:,:),allocatable, save :: Stor3Up
  double precision, dimension(:,:),allocatable, save :: Stor4Up
  double precision, dimension(:,:),allocatable, save :: Stor1Lo
  double precision, dimension(:,:),allocatable, save :: Stor2Lo
  double precision, dimension(:,:),allocatable, save :: Stor3Lo
  double precision, dimension(:,:),allocatable, save :: Stor4Lo
  double precision, dimension(:,:),allocatable, save :: Stor1m
  double precision, dimension(:,:),allocatable, save :: Stor2m
  double precision, dimension(:,:),allocatable, save :: Stor3m
  double precision, dimension(:,:),allocatable, save :: Stor4m
  double precision, dimension(:,:),allocatable, save :: Stor1z
  double precision, dimension(:,:),allocatable, save :: Stor2z
  double precision, dimension(:,:),allocatable, save :: Stor3z
  double precision, dimension(:,:),allocatable, save :: Stor4z
  double precision, dimension(:,:),allocatable, save :: PP1
  double precision, dimension(:,:),allocatable, save :: PP2
  double precision, dimension(:,:),allocatable, save :: DJ1
  double precision, dimension(:,:),allocatable, save :: DJ2
  double precision, dimension(:,:),allocatable, save :: FD1
  double precision, dimension(:,:),allocatable, save :: FD2
  double precision, dimension(:,:),allocatable, save :: DJP1
  double precision, dimension(:,:),allocatable, save :: DJP2
  double precision, dimension(:,:),allocatable, save :: DJP3
  double precision, dimension(:,:),allocatable, save :: DJP4
  !
  !
  !      Variables that don't need to be stored for subsequent
  !      right-hand sides.
  !
  double precision f1(1:m1,1:m2)
  double precision f2(1:m1,1:m2),f3(1:m2,1:m1),f4(1:m2,1:m2)
  double precision A1(1:m1,1:3),B1(1:m1,1:5),C1(1:m1,1:3)
  double precision A2(1:m2,1:3),B2(1:m2,1:5),C2(1:m2,1:3)
  double precision Qa1(1:m1,1:2),Qa2(1:m2,1:2)
  double precision W1(1:2,1:m1),W2(1:2,1:m2)
  double precision konst1(1:m1),konst2(1:m2)
  double precision Tut1(1:m1),Tut2(1:m2)
  double precision DB1(1:m1),DB2(1:m2)
  double precision sh(0:N-5,0:N-5)
  double precision norm1,norm2,norm3,norm4
  !
  !      Variable for den iterative delen
  !
  double precision ka1(1:m1),ka2(1:m1),ka3(1:m1),ka4(1:m1)
  double precision la1(1:m2),la2(1:m2),la3(1:m2),la4(1:m2)
  !
  call dtime(tt)
!  count = 0
!  count_max = 0
!  count_rate = 0
!  call system_clock(count,count_rate,count_max)
!  ocount = count

  !      flag = 0     vil si at en har foerste hoyreside, ellers
  !                   er det en paafoelgende hoyreside.

  if (flag .EQ. 0 .AND. oldn .NE. 0) then
     !
     !     We have been called before, but with a different N, must
     !     deallocate all arrays before allocating new ones.
     !
     deallocate (rr1)
     deallocate (rr2)
     deallocate (rr3)
     deallocate (ss2)
     deallocate (ss3)
     deallocate (D1)
     deallocate (p1)
     deallocate (q1)
     deallocate (r1)
     deallocate (s1)
     deallocate (D2)
     deallocate (p2)
     deallocate (q2)
     deallocate (r2)
     deallocate (s2)
     deallocate (dd)
     deallocate (e1)
     deallocate (e2)
     deallocate (e3)
     deallocate (e4)
     deallocate (Id1)
     deallocate (Id1inv)
     deallocate (Id2)
     deallocate (Id2inv)
     deallocate (Ja1)
     deallocate (Fa1)
     deallocate (Ja2)
     deallocate (Fa2)
     deallocate (wsave)
     deallocate (Sum1)
     deallocate (Sum2)
     deallocate (Stor1Up)
     deallocate (Stor2Up)
     deallocate (Stor3Up)
     deallocate (Stor4Up)
     deallocate (Stor1Lo)
     deallocate (Stor2Lo)
     deallocate (Stor3Lo)
     deallocate (Stor4Lo)
     deallocate (Stor1m)
     deallocate (Stor2m)
     deallocate (Stor3m)
     deallocate (Stor4m)
     deallocate (Stor1z)
     deallocate (Stor2z)
     deallocate (Stor3z)
     deallocate (Stor4z)
     deallocate (Sum1e1)
     deallocate (Sum2e1)
     deallocate (Sum3e1)
     deallocate (Sum4e1)
     deallocate (Sum1e2)
     deallocate (Sum2e2)
     deallocate (Sum3e2)
     deallocate (Sum4e2)
     deallocate (Sum1e3)
     deallocate (Sum2e3)
     deallocate (Sum3e3)
     deallocate (Sum4e3)
     deallocate (Sum1e4)
     deallocate (Sum2e4)
     deallocate (Sum3e4)
     deallocate (Sum4e4)
     deallocate (Diag1)
     deallocate (Diag2)
     deallocate (Dinv1)
     deallocate (Dinv2)
     deallocate (PP1)
     deallocate (PP2)
     deallocate (DJ1)
     deallocate (DJ2)
     deallocate (FD1)
     deallocate (FD2)
     deallocate (DJP1)
     deallocate (DJP2)
     deallocate (DJP3)
     deallocate (DJP4)
     !
  end if
  !
  !  dealloate done
  !
  oldn = n
  !
  if (flag .EQ. 0) then
     !
     !
     !      Allocate the necessary variables, only to be done for the
     !      first right hand side.
     !
     !
     allocate (rr1(0:N-5))
     allocate (rr2(0:N-5))
     allocate (rr3(0:N-5))
     allocate (ss2(0:N-5))
     allocate (ss3(0:N-5))
     allocate (D1(1:m1))
     allocate (p1(1:m1))
     allocate (q1(1:m1))
     allocate (r1(1:m1))
     allocate (s1(1:m1))
     allocate (D2(1:m2))
     allocate (p2(1:m2))
     allocate (q2(1:m2))
     allocate (r2(1:m2))
     allocate (s2(1:m2))
     allocate (dd(0:N-5))
     allocate (e1(1:m1**2))
     allocate (e2(1:m1*m2))
     allocate (e3(1:m1*m2))
     allocate (e4(1:m2**2))
     !
     allocate (Id1(1:m1,1:m1))
     allocate (Id1inv(1:m1,1:m1))
     allocate (Id2(1:m2,1:m2))
     allocate (Id2inv(1:m2,1:m2))
     allocate (Ja1(1:m1,1:2))
     allocate (Fa1(1:2,1:m1))
     allocate (Ja2(1:m2,1:2))
     allocate (Fa2(1:2,1:m2))
     allocate (Sum1(1:m1,1:m1))
     allocate (Sum2(1:m2,1:m2))
     allocate (Stor1Up(1:m1,1:3*m1))
     allocate (Stor2Up(1:m2,1:3*m1))
     allocate (Stor3Up(1:m1,1:3*m2))
     allocate (Stor4Up(1:m2,1:3*m2))
     allocate (Stor1Lo(1:m1,1:3*m1))
     allocate (Stor2Lo(1:m2,1:3*m1))
     allocate (Stor3Lo(1:m1,1:3*m2))
     allocate (Stor4Lo(1:m2,1:3*m2))
     allocate (Stor1m(1:m1,1:m1))
     allocate (Stor2m(1:m2,1:m1))
     allocate (Stor3m(1:m1,1:m2))
     allocate (Stor4m(1:m2,1:m2))
     allocate (Stor1z(1:m1,1:m1))
     allocate (Stor2z(1:m2,1:m1))
     allocate (Stor3z(1:m1,1:m2))
     allocate (Stor4z(1:m2,1:m2))
     allocate (Sum1e1(1:m1))
     allocate (Sum2e1(1:m1))
     allocate (Sum3e1(1:m1))
     allocate (Sum4e1(1:m1))
     allocate (Sum1e2(1:m2))
     allocate (Sum2e2(1:m2))
     allocate (Sum3e2(1:m2))
     allocate (Sum4e2(1:m2))
     allocate (Sum1e3(1:m1))
     allocate (Sum2e3(1:m1))
     allocate (Sum3e3(1:m1))
     allocate (Sum4e3(1:m1))
     allocate (Sum1e4(1:m2))
     allocate (Sum2e4(1:m2))
     allocate (Sum3e4(1:m2))
     allocate (Sum4e4(1:m2))
     allocate (Diag1(1:m1))
     allocate (Diag2(1:m2))
     allocate (wsave(1:6*N+27))
     allocate (Dinv1(1:m1))
     allocate (Dinv2(1:m2))
     allocate (PP1(1:2,1:2))
     allocate (PP2(1:2,1:2))
     allocate (DJ1(1:m1,1:2))
     allocate (DJ2(1:m2,1:2))
     allocate (FD1(1:2,1:m1))
     allocate (FD2(1:2,1:m2))
     allocate (DJP1(1:2*m1,1:4))
     allocate (DJP2(1:2*m2,1:4))
     allocate (DJP3(1:2*m1,1:4))
     allocate (DJP4(1:2*m2,1:4))
     !
     !
     !      The spectral matrices A1,A2,B1,B2,C1,C2 are constructed.
     !      Only the 3 main diagonals of the upper triangular matrices
     !      A1 and A2 are computed, the rest of the elements of A_i can
     !      be computed from p_i,q_i,r_i and s_i for i = 1,2.
     !
     !      A_i^(-1)B_i-(A_i^(-1)C_i)^2=Qa_i*W_i for i = 1,2
     !
     !
     odde = 1
     call cheb_disk(m1,A1,B1,C1,D1,p1,q1,r1,s1,Qa1,W1,odde)
     odde = 0
     call cheb_disk(m2,A2,B2,C2,D2,p2,q2,r2,s2,Qa2,W2,odde)  
     !
     !
     !     The eigenvalue problems Idinv_i*A_i^(-1)*Id_i = Diag_i
     !     are solved in O(N^2) flops.
     !
     !
     odde = 1
     call cheb_eig(m1,A1,D1,p1,q1,r1,s1,C1,odde,Id1,Id1inv,Diag1)
     odde = 0
     call cheb_eig(m2,A2,D2,p2,q2,r2,s2,C2,odde,Id2,Id2inv,Diag2)
     !
     !
     !     Ja_i = Idinv_i*Qa_i and Fa_i = W_i*Id_i are computed.
     !
     !
     do i = 1,m1
        Ja1(i,1) = Id1inv(i,1)*Qa1(1,1)
        Ja1(i,2) = dot_product(Id1inv(i,:),Qa1(:,2))
        Fa1(1,i) = dot_product(W1(1,:),Id1(:,i))
        Fa1(2,i) = W1(2,m1)*Id1(m1,i)
     end do
     !
     do i = 1,m2
        Ja2(i,1) = Id2inv(i,1)*Qa2(1,1)
        Ja2(i,2) = dot_product(Id2inv(i,:),Qa2(:,2))
        Fa2(1,i) = dot_product(W2(1,:),Id2(:,i))
        Fa2(2,i) = W2(2,m2)*Id2(m2,i)
     end do
     !
     DB1 = Diag1**2
     DB2 = Diag2**2
     !
     !
     !     Factorization of the block diagonal matrix hatZ from sec. 4.6.2.
     !
     !
     call cheb_sepfakt(A1,B1,C1,DB1,Diag1,p1,q1,r1,s1,Stor1Lo,     &
          Stor1Up,Stor1z,Stor1m,alf,bet,m1,m1)
     !
     call cheb_sepfakt(A2,B2,C2,DB1,Diag1,p2,q2,r2,s2,Stor2Lo,     &
          Stor2Up,Stor2z,Stor2m,alf,bet,m2,m1)
     !
     call cheb_sepfakt(A1,B1,C1,DB2,Diag2,p1,q1,r1,s1,Stor3Lo,     &
          Stor3Up,Stor3z,Stor3m,alf,bet,m1,m2)
     !
     call cheb_sepfakt(A2,B2,C2,DB2,Diag2,p2,q2,r2,s2,Stor4Lo,     &
          Stor4Up,Stor4z,Stor4m,alf,bet,m2,m2)
     !
     Tut1 = 0.5d0*(1.d0+alf*Diag1)
     Tut2 = 0.5d0*(1.d0+alf*Diag2)
     !
     konst1 = 0.5d0*(1.d0+bet*DB1+alf*Diag1)
     konst2 = 0.5d0*(1.d0+bet*DB2+alf*Diag2)
     !
     !
     !     Construction of the matrices D^(k) from sec. 4.6.3.
     !
     !
     call compDiag(konst1,Diag1,Diag1,DB1,DB1,Tut1,e1,m1,m1)
     call compDiag(konst1,Diag2,Diag1,DB2,DB1,Tut2,e2,m2,m1)
     call compDiag(konst2,Diag1,Diag2,DB1,DB2,Tut1,e3,m1,m2)
     call compDiag(konst2,Diag2,Diag2,DB2,DB2,Tut2,e4,m2,m2)
     !
     !
     !     The sums that will be used in the iteretive solution of P
     !     are computed.
     !
     !
     ka1 = Ja1(:,1)*Fa1(1,:)
     ka2 = Ja1(:,2)*Fa1(1,:)
     ka3 = Ja1(:,1)*Fa1(2,:)
     ka4 = Ja1(:,2)*Fa1(2,:)
     !
     la1 = Ja2(:,1)*Fa2(1,:)
     la2 = Ja2(:,2)*Fa2(1,:)
     la3 = Ja2(:,1)*Fa2(2,:)
     la4 = Ja2(:,2)*Fa2(2,:)
     !
     call cheb_compSum(ka1,ka2,ka3,ka4,e1,Sum1e1,Sum2e1,Sum3e1,Sum4e1,m1,m1)
     call cheb_compSum(ka1,ka2,ka3,ka4,e2,Sum1e2,Sum2e2,Sum3e2,Sum4e2,m2,m1)
     call cheb_compSum(la1,la2,la3,la4,e3,Sum1e3,Sum2e3,Sum3e3,Sum4e3,m1,m2)
     call cheb_compSum(la1,la2,la3,la4,e4,Sum1e4,Sum2e4,Sum3e4,Sum4e4,m2,m2)
     !
     !
     !     Also matrices needed for the preconditioning, Dinv, DJ, FD
     !     PP are computed.
     !
     !
     call cg_cheb_pre_pfirst(alf,bet,Ja1,Fa1,Diag1,Dinv1,DJ1,FD1, &
          PP1,m1)
     !
     call cg_cheb_pre_pfirst(alf,bet,Ja2,Fa2,Diag2,Dinv2,DJ2,FD2, &
          PP2,m2)
     !
     !     Compute the first part of P, E_i^(-1) G_i^(-1) E_i, i=1,2.
     !     Need not be called for bet=0, (see 4.70)
     !
     !
     if ((bet .gt. 0.d0) .or. (bet .lt. 0.d0)) then
        !
        call cg_cheb_pfirst(Sum1e1,Sum2e1,    &
             Sum3e1,Sum4e1,Dinv1,DJ1,FD1,PP1,DJP1,m1)

        call cg_cheb_pfirst(Sum1e2,Sum2e2,    &
             Sum3e2,Sum4e2,Dinv2,DJ2,FD2,PP2,DJP2,m2)

        call cg_cheb_pfirst(Sum1e3,Sum2e3,    &
             Sum3e3,Sum4e3,Dinv1,DJ1,FD1,PP1,DJP3,m1)

        call cg_cheb_pfirst(Sum1e4,Sum2e4,    &
             Sum3e4,Sum4e4,Dinv2,DJ2,FD2,PP2,DJP4,m2)
        !
     end if
     !
     !
     !     Initialization of vr2cos and vfrcos.
     !
     !
     call vinifft(N-1,wsave)
     !
     !
     !     The coefficients needed for going from Shen's basis to
     !     standard Chebyshev basis are computed.
     !
     !
     do i = 0,N-5
        rk = i
        dd(i) = 1.d0
        ss2(i) = -2*(rk+2)/(rk+3)
        ss3(i) = (rk+1)/(rk+3)
        rr1(i) = pi/2
        rr2(i) = ss2(i)*pi/2
        rr3(i) = ss3(i)*pi/2
     end do
     !
     rr1(0)=2*rr1(0)
     !
  endif
  !
  !
  !     End of pre-prosessing for a first right hand side.
  !
  !
  !     From here for repeated solves.
  !
  !
  !     Preprocessing time
  !
  call dtime(tt)
  print*,'UserTime for preprocessing: ',tt(1)
!  call system_clock(count,count_rate,count_max)
!  tid1 = (count-ocount)/(1.d0 * count_rate)
!  ocount = count
  call dtime(tt)
  !
  !
  !     Transformation from physical space to spectral space, using
  !     FFT's.
  !
  !      
  flagget = 0
  call vr2cos(hoyre,N-1,flagget,wsave)
  !
!  call system_clock(count,count_rate,count_max)
!  tid3=(count-ocount)/(1.d0 * count_rate)
!  ocount = count
  !
  !
  !     Transformation from standard Chebyshev basis to Shen's special
  !     basis.
  !
  !
  hoyre(0:N-5,0:N-5) = leg_to_shen(hoyre,dd,rr1,rr2,rr3,N)
  !
  !
  !     The right-hand side is split into four parts.
  !
  !
    f1 = matr_to_vekt(hoyre(0:N-5,0:N-5),m1,m1,1,1,N)
    f2 = matr_to_vekt(hoyre(0:N-5,0:N-5),m2,m1,0,1,N)
    f3 = matr_to_vekt(hoyre(0:N-5,0:N-5),m1,m2,1,0,N)
    f4 = matr_to_vekt(hoyre(0:N-5,0:N-5),m2,m2,0,0,N)

!  call matr_to_sub(hoyre(0:N-5,0:N-5),f1,m1,m1,1,1,N)
!  call matr_to_sub(hoyre(0:N-5,0:N-5),f2,m2,m1,0,1,N)
!  call matr_to_sub(hoyre(0:N-5,0:N-5),f3,m1,m2,1,0,N)
!  call matr_to_sub(hoyre(0:N-5,0:N-5),f4,m2,m2,0,0,N)
  !
  !
  !     f_i = kron(A^(-1),I)*f_i  i = 1,..,4
  !
  !
  call  Amvprod(D1,p1,q1,r1,s1,f1,m1,m1)
  call  Amvprod(D1,p1,q1,r1,s1,f2,m2,m1)
  call  Amvprod(D2,p2,q2,r2,s2,f3,m1,m2)
  call  Amvprod(D2,p2,q2,r2,s2,f4,m2,m2)
  !
  !
  !    Norms of f1,...,f4 are computed to test for zero right hand side.
  !
  !
  !  norm1 = sqrt(dot_product(f1,f1))
  !  norm2 = sqrt(dot_product(f2,f2))
  !  norm3 = sqrt(dot_product(f3,f3))
  !  norm4 = sqrt(dot_product(f4,f4))
  !
  !
  !     The four systems M_1 u_1 = f_1 , ..., M_4 u_4 = f_4 from
  !     (4.56)-(4.59) are solved using cg_cheb_hele.
  !
  !
  !  if (norm1.gt.rhs_norm) then
  call cg_cheb_hele(m1,m1,Diag1,Id1,Id1inv,Id1,Id1inv,D1,p1,q1,r1,s1,  &
       f1,Ja1,Ja1,Fa1,Fa1,e1, &
       Sum1e1,Sum2e1,Sum3e1,Sum4e1,Stor1Up,Stor1Lo,Stor1z,Stor1m,  &
       Dinv1,PP1,DJ1,FD1,DJP1,alf,bet)
  ! else
  !     f1 = 0.d0
  !  end if
  !
  !
  !  if (norm2.gt.rhs_norm) then
  call cg_cheb_hele(m2,m1,Diag1,Id1,Id1inv,Id2,Id2inv,D2,p2,q2,r2,s2,   &
       f2,Ja2,Ja1,Fa2,Fa1,e2, &
       Sum1e2,Sum2e2,Sum3e2,Sum4e2,Stor2Up,Stor2Lo,Stor2z,Stor2m,  &
       Dinv2,PP2,DJ2,FD2,DJP2,alf,bet)
  !  else
  !     f2 = 0.d0
  !  end if
  !
  !  if (norm3.gt.rhs_norm) then
  call cg_cheb_hele(m1,m2,Diag2,Id2,Id2inv,Id1,Id1inv,D1,p1,q1,r1,s1,  &
       f3,Ja1,Ja2,Fa1,Fa2,e3, &
       Sum1e3,Sum2e3,Sum3e3,Sum4e3,Stor3Up,Stor3Lo,Stor3z,Stor3m,  &
       Dinv1,PP1,DJ1,FD1,DJP3,alf,bet)
  !  else
  !     f3 = 0.d0
  !  end if
  !
  !  if (norm4.gt.rhs_norm) then
  call  cg_cheb_hele(m2,m2,Diag2,Id2,Id2inv,Id2,Id2inv,D2,p2,q2,r2,s2, &
       f4,Ja2,Ja2,Fa2,Fa2,e4, &
       Sum1e4,Sum2e4,Sum3e4,Sum4e4,Stor4Up,Stor4Lo,Stor4z,Stor4m,  &
       Dinv2,PP2,DJ2,FD2,DJP4,alf,bet)
  !  else
  !     f4 = 0.d0
  !  end if
  !
  !
  !     The four solutions are merged together.
  !
  !
    call vekt_to_matr(f1,sh,m1,m1,1,1,N)
    call vekt_to_matr(f2,sh,m2,m1,0,1,N)
    call vekt_to_matr(f3,sh,m1,m2,1,0,N)
    call vekt_to_matr(f4,sh,m2,m2,0,0,N)

  !call sub_to_matr(f1,sh,m1,m1,1,1,N)
  !call sub_to_matr(f2,sh,m2,m1,0,1,N)
  !call sub_to_matr(f3,sh,m1,m2,1,0,N)
  !call sub_to_matr(f4,sh,m2,m2,0,0,N)
  !
  !
  !     The solution is transformed from Shen's basis back to
  !     standard Chebyshev basis.
  !
  !
  call basistrans(hoyre,sh,dd,ss2,ss3,N)
  !
!  call system_clock(count,count_rate,count_max)
!  tid2 = (count-ocount)/(1.d0 * count_rate)
!  ocount = count
  !
  !
  !     Transform from spectral to physical space
  !
  !
  flagget = 1
  call vr2cos(hoyre,N-1,flagget,wsave)
  call dtime(tt)
  print*,'UserTime for Transform and Solution: ',tt(1)
!  call system_clock(count,count_rate,count_max)
!  tid3 = tid3 + (count-ocount)/(1.d0 * count_rate)
  !
!  tid4 = tid1+tid2+tid3
  !
end subroutine cg_cheb

