1C> \ingroup selci
2C> @{
3C>
4      subroutine selci_eijil(e,i,j,l,ns,indbar,iocc,
5     $     w1,w2,work,numf,numf2)
6*
7* $Id$
8*
9#include "implicit.fh"
10#include "errquit.fh"
11#include "ciinfo.fh"
12      dimension e(*),indbar(*),iocc(*),work(*),
13     $     w1(nfmax,nfmax,nsmax),w2(nfmax2,nfmax,nsmax-1)
14c
15c     e(u,v) = <Iu|Eij,il|Iv>, where u,v label the spin functions
16c     and the orbital occupation of I is specified by ns, iocc
17c     and indbar. Only for i.ne.j.ne.l. j and l are assumed
18c     to be in order of increasing occupation in I
19c
20c     Note that this value is independent of i and can be simplified later
21c
22c     possible cases are:
23c
24c     iocc(j)  iocc(l)  case
25c       0         0       1
26c       0         1       2
27c       1         1       3
28c
29      if (i.eq.j .or. i.eq.l .or. j.eq.l)
30     $     call errquit('eijil: one of i=j,i=l,j=l',i, UNKNOWN_ERR)
31      if (iocc(j).gt.iocc(l))
32     $     call errquit('eijil: iocc ordering',iocc(j), UNKNOWN_ERR)
33      numf = nf(ns)
34      numf2 = nf(ns+2)
35      ib = indbar(i) + (nsmax-ns-2)
36*     call selci_mxma(w2(1,1,ib),1,nfmax2,w1(1,1,ib),1,nfmax,
37*     $     e,1,numf,numf,numf2,numf2)
38      call selci_axb(w2(1,1,ib),nfmax2,w1(1,1,ib),nfmax,
39     $     e,numf,numf,numf2,numf2)
40c
41      icase = iocc(j) + iocc(l) + 1
42      goto (10,20,30) icase
43      call errquit('eijil: invalid case',icase, UNKNOWN_ERR)
44c
45c     nj=0, nl=0: e(u,v) = eiaeib(u,p)*w1(q,p,lb)*w1(v,q,jb)
46c                                    ab     ->  la   ->     jl
47c
48 10   continue
49      lb = indbar(l) + (nsmax-ns-2)
50      jb = indbar(j) + (nsmax-ns-2)
51      if (j.gt.l) jb = jb + 1
52*      call selci_mxma(e,1,numf,w1(1,1,lb),nfmax,1,
53*     $     work,1,numf,numf,numf2,numf2)
54      call selci_axbt(e,numf,w1(1,1,lb),nfmax,
55     $     work,numf,numf,numf2,numf2)
56*      call selci_mxma(work,1,numf,w1(1,1,jb),nfmax,1,
57*     $     e,1,numf,numf,numf2,numf2)
58      call selci_axbt(work,numf,w1(1,1,jb),nfmax,
59     $     e,numf,numf,numf2,numf2)
60      return
61c
62c     nj=0, nl=1: e(u,v) = eiaeib(u,p)*w2(q,p,lb)*w1(v,q,jb)
63c                                   lab     -> lla   ->    llj
64c
65 20   continue
66      lb = indbar(l) + (nsmax-ns-2)
67      jb = indbar(j) + (nsmax-ns)
68      if (j.gt.l) jb = jb - 1
69*      call selci_mxma(e,1,numf,w2(1,1,lb),nfmax2,1,
70*     $     work,1,numf,numf,numf2,numf)
71      call selci_axbt(e,numf,w2(1,1,lb),nfmax2,
72     $     work,numf,numf,numf2,numf)
73*      call selci_mxma(work,1,numf,w1(1,1,jb),nfmax,1,
74*     $     e,1,numf,numf,numf,numf)
75      call selci_axbt(work,numf,w1(1,1,jb),nfmax,
76     $     e,numf,numf,numf,numf)
77      numf2 = numf
78      return
79c
80c     nj=1, nl=1: e(u,v) = eiaeib(u,p)*w2(q,p,lb)*w1(v,q,jb)
81c                                  jlab    -> llja   ->   jjll
82c
83 30   continue
84      lb = indbar(l) + (nsmax-ns-2)
85      jb = indbar(j) + (nsmax-ns)
86      if (j.gt.l) jb = jb - 1
87      numf2m = nf(ns-2)
88*      call selci_mxma(e,1,numf,w2(1,1,lb),nfmax2,1,
89*     $     work,1,numf,numf,numf2,numf)
90      call selci_axbt(e,numf,w2(1,1,lb),nfmax2,
91     $     work,numf,numf,numf2,numf)
92*      call selci_mxma(work,1,numf,w2(1,1,jb),nfmax2,1,
93*     $     e,1,numf,numf,numf,numf2m)
94      call selci_axbt(work,numf,w2(1,1,jb),nfmax2,
95     $     e,numf,numf,numf,numf2m)
96      numf2 = numf2m
97c
98      end
99C>
100C> @}
101