
c/* %W%      latest revision %G% %U% */

      subroutine Fmatrx( den , n1 , nmax , nspin , nptmax , ldum, psfac)
********************************************************************************

      implicit real*8 ( a-h, o-z )
      complex*16 zf, zwork, zi
      integer ipiv, lwork, info, lda

c     dimension   ur(97,97,8) , ui(97,97,8) , zf(194,194), zwork(100)
      dimension   ur(97,97,8) , ui(97,97,8) , zf(86,86), zwork(100)
      dimension   nifty(20), ipiv(106)   
      dimension   den(*)

      common /nlspfl/   ur, ui
      common /foptp/    zf
      common /switch/   nifty
      common /spins/    nspina , nspinb , nspinc , nspind , nspine

********************************************************************************

         zi = (0.,1.0)
         lwork = 100
c        write(11,666) ldum, nspin, ur(n1,n1,nspin), ui(n1,n1,nspin)
 666     format(2i4, 2e14.7)
       si = psfac
*******************************************************************
      do 240 i2 = 1,nmax
         do 240 i3 = 1,nmax

c ***
c *** spin 1/2 x 1/2     do nspin=3/4  5/6 together
c ***
c ***       |ra|   |va|     |va vd|   |ra|
c ***       |  | = |  |  +  |     | * |  |
c ***       |rb|   |vb|     |vb vc|   |rb|
c ***
c *** (abcd) = (3456)for j=l+1(nspin=4)
c ***        = (5634)for j=l-1 (nspin=6)

 221        if (i2 .le. n1) go to 224
            ibc = i2 - n1
            if (i3 .le. n1) go to 222

c *** element c

            jc = i3 - n1
c changed index of den to jc 15 Dec 88
            zf(i2,i3) = ur(ibc,jc,nspinc) * den(jc)
     &                + zi * ui(ibc,jc,nspinc) * den(jc) 
********************* New T-Mat 7/27/92  tvm ********************************
            if (jc.eq.n1) then
            zf(i2,i3) = zf(i2,i3) - si * ui(ibc,jc,nspinc)
     &                  + si * zi * ur(ibc,jc,nspinc)
            endif
*****************************************************************************

            go to 229
  224       if (i3 .le. n1) go to 223

c *** element d

            jd = i3 - n1
c changed index of den to jd 15 Dec 88
            zf(i2,i3) = ur(i2,jd,nspind) * den(jd)
     &                 + zi * ui(i2,jd,nspind) * den(jd)
********************* New T-Mat 7/27/92  tvm ********************************
            if (jd.eq.n1) then
            zf(i2,i3) = zf(i2,i3) - si * ui(i2,jd,nspind)
     &                   + zi * si * ur(i2,jd,nspind)
            endif
*****************************************************************************

            go to 229

c *** element b

  222       zf(i2,i3) = ur(ibc,i3,nspinb) * den(i3)
     &                   + zi * ui(ibc,i3,nspinb) * den(i3)
********************* New T-Mat 7/27/92  tvm ********************************
            if (i3.eq.n1) then
            zf(i2,i3) = zf(i2,i3) - si * ui(ibc,i3,nspinb)
     &                  + zi * si * ur(ibc,i3,nspinb)
            endif
*****************************************************************************

            go to 229

c *** element a

 223        zf(i2,i3) = ur(i2,i3,nspina) * den(i3)
     &                   + zi * ui(i2,i3,nspina) * den(i3)
********************* New T-Mat 7/27/92  tvm ********************************
            if (i3.eq.n1) then
            zf(i2,i3) = zf(i2,i3) - si * ui(i2,i3,nspina)
     &                  + zi * si * ur(i2,i3,nspina)
            endif
*****************************************************************************

 229        if (i2 .ne. i3) go to 240

c *** diagonal term

            zf(i2,i3) = zf(i2,i3) + 1.
 240  continue

c *** cmatin  does  matrix-inversion.

c *** write out at least part of f matrix

      nfdim = nptmax * 2

c *** cmatin does matrix inversion for copmplex matrix

c     call CMATIN( fr , fi , nmax , nfdim )
      call zgetrf( nmax, nmax, zf, nmax, ipiv, info)
      call zgetri( nmax, zf, nmax, ipiv, zwork, lwork, info)
      if (info.ne.0) then
      write(6,2323) info
      endif
 2323  format(22hzgetri didn't work!!! ,i6)
c     write(7,*) zwork(1)
c     write(7,1492) info
 1492 format(i9)

      RETURN

 250  nms = 1 + nifty(6)
      write(6,820) nms

      RETURN

c *** Formats ***

 227  format(i3, 8e15.6 / 200(3x, 8e15.6 / ))
 820  format (' R matrix calculated with ms series up to order ',i3)
 928  format('0',' re/im f matrix before inversion ')
 1498 format('  ???fmatrx  ur(lid,lid,1)= ',e10.4,'   ui(=)= ',e10.4)

      end
