c23456789012345678901234567890123456789012345678901234567890123456789012


      SUBROUTINE alphaxx(kcomp) 

c     ******************************************************************
c
c     This routine computes the vector alpha 
c
c     ******************************************************************

c     -------------
c     Declarations
c     -------------

c     ''''''''''''''''''''''
      implicit none
      include 'commonvar.f'      
c     ''''''''''''''''''''''

c---- Argument variables

      integer	kcomp


c---- Local variables

      integer	n,n1,n2,n3,kco

      real*8	MAT_alpha(idim_n_nl)	

      real*8	FFT(idim_n_o,idim_n_nl)

      character date*8,time*10,zone*5

      integer   values(8),ipvt(idim_n_nl)




c     --------------------
c     Matrix F Coa F' + Cf
c     --------------------

c---- Product of F and Coa 

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) '  '//time(1:2)//':'//
     >                 time(3:4)//':'//
     >                 time(5:6)//'  Get F x Coa ...'

      MAT2(1:n_nl,1:n_o) = matmul(FF (1:n_nl,1:n_o),
     >				  Coa(1:n_o ,1:n_o))


c---- Product of (F Coa) and beta

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) '  '//time(1:2)//':'//
     >                 time(3:4)//':'//
     >                 time(5:6)//'  Get F x Coa x beta ...'

      vec1(1:n_nl) = matmul(MAT2(1:n_nl,1:n_o),
     >		            beta(       1:n_o))


c---- Product of (F Coa) and F' 

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) '  '//time(1:2)//':'//
     >                 time(3:4)//':'//
     >                 time(5:6)//'  Get F x Coa x FT ...'

      do n1=1,n_o
      do n2=1,n_nl
         FFT(n1,n2) = FF(n2,n1)
      enddo
      enddo

      MAT1(1:n_nl,1:n_nl) = matmul(MAT2(1:n_nl,1:n_o),
     >			           FFT (1:n_o ,1:n_nl))


c---- Matrix F Coa F'  +   Cf

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) '  '//time(1:2)//':'//
     >                 time(3:4)//':'//
     >                 time(5:6)//'  Get F Coa FT   +   Cf ...' 

      do n=1,n_nl
         MAT1(n,n) = MAT1(n,n) + 1./Cf_1(n)
      enddo




c     ----------------------------------------------------------
c     Solve the linear system (F Coa FT + Cf) alpha = F Coa beta 
c
c     Note
c     ----
c     LU decomposition is used
c     ----------------------------------------------------------

c---- Save MAT1 and vec1 

      do n1=1,n_nl
      do n2=1,n_nl
         MAT2(n1,n2) = MAT1(n1,n2)
      enddo
      enddo
c--
      do n=1,n_nl
         vec2(n) = vec1(n)
      enddo


c---- LU factorize F Coa FT + Cf = L x U ...' using ESSL routine

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) '  '//time(1:2)//':'//
     >                 time(3:4)//':'//
     >                 time(5:6)//'  LU factorize F Coa FT + Cf ...'

c     ----------------------------------
      call DGEF(MAT1,idim_n_o,n_nl,ipvt)	 
c     ----------------------------------


c---- Solve (F Coa FT + Cf) alpha   =   F Coa beta' using ESSL routine

c     -----------------------------------------
      call date_and_time(date,time,zone,values)
c     -----------------------------------------
      write(*,*) '  '//time(1:2)//':'//
     >                 time(3:4)//':'//
     >                 time(5:6)//'  Solve LU alpha = F Coa beta ...'

c     -----------------------------------------
      call DGES(MAT1,idim_n_o,n_nl,ipvt,vec1,0)        
c     -----------------------------------------




c     --------------------------------
c     Assignation of the vector alpha
c     --------------------------------

      do n=1,n_nl
         alpha(n) = vec1(n)
      enddo




c     -----------
c     Diagnostics 
c     -----------

c---- RETURN statement

      return


c---- Compute MAT1 x alpha 

      write(*,*) '  Computing solution diagnostics (in check2) ...'

      do n1=1,n_nl
         MAT_alpha(n1) = 0.d0
         do n2=1,n_nl
            MAT_alpha(n1) =
     >      MAT_alpha(n1) + MAT2(n1,n2)*alpha(n2)
         enddo
      enddo


c---- Writing diagnostics

      if (kcomp.eq.1) open (50,file='check2.'//tag)
      do n=1,n_nl
         write(50,*) n,MAT_alpha(n)-
     >                 vec2     (n)
      enddo




c     -------------------
c     End of the routine
c     -------------------

      return

      end

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

