1C$Id$
2************************************************************************
3*                                                                      *
4      subroutine ecp_cstrans (l,ncart,ntr,l_hi,l_lo,nsph,cart,ldc,
5     &    sph,lds,csco,lcsco,lstart,indicator,position,mode)
6*                                                                      *
7*   Transform from cartesian monomials x^i y^j z^k to real solid       *
8*   spherical tensor functions r^n G_lm(theta,phi). These are          *
9*   sqrt(2pi(1+delta_m,0)/2l+1) times the real spherical harmonics.    *
10*                                                                      *
11*   l (inp) - angular momentum of shell, l = i+j+k                     *
12*   ncart (inp) - number of cartesian components = (l+1)*(l+2)/2       *
13*   ntr (inp) - number of quantities to transform                      *
14*   l_hi - maximum angular momentum of spherical quantities            *
15*   l_lo - minimum angular momentum of spherical quantities            *
16*   nsph (out) - number of spherical components                        *
17*   cart (i/o) - array of cartesian quantities.                        *
18*   ldc (inp) - leading dimension of cart                              *
19*   sph (i/o) - array of spherical quantities                          *
20*   lds (inp) - leading dimension of sph                               *
21*   csco (inp) - transformation arrays                                 *
22*   lcsco (inp) - length of transformation arrays                      *
23*   lstart - pointers to beginning of transformation for given l       *
24*   indicator (inp) - +1 = cartesian to spherical,                     *
25*                     -1 = spherical to cartesian                      *
26*   position (inp) - +1 = post-multiply, array is ntr*(ncart or nsph)  *
27*                    -1 = premultiply, array is (ncart or nsph)*ntr    *
28*   mode (inp) - +1 = transform with regular matrix,                   *
29*                -1 = transform with inverse matrix                    *
30*                                                                      *
31*   Written by K. G. Dyall                                             *
32*                                                                      *
33************************************************************************
34      implicit none
35#include "errquit.fh"
36      integer l,ncart,ntr,l_hi,l_lo,nsph,ldc,lds,lcsco,lstart(0:1,0:l),
37     &    indicator,position,mode
38      integer i,j,offset_start,offset_end
39      double precision cart(ldc,ncart),sph(lds,ncart)
40      double precision csco(lcsco),zero,one
41      parameter (zero = 0.0d00, one = 1.0d00)
42*
43*   Test validity of input parameters
44*
45      if (abs(indicator) .ne. 1) call errquit (
46     &      'Illegal value of indicator in ecp_cstrans',99, BASIS_ERR)
47      if (abs(position) .ne. 1) call errquit (
48     &      'Illegal value of position in ecp_cstrans',99, BASIS_ERR)
49      if (abs(mode) .ne. 1) call errquit (
50     &      'Illegal value of mode in ecp_cstrans',99, BASIS_ERR)
51*
52C      write (6,*) 'l_hi,l_lo',l_hi,l_lo
53      offset_start = (l_hi+1)*(l_hi+2)/2
54      offset_end = (l_lo-1)*(l_lo)/2
55C      write (6,*) 'offset_start,offset_end',offset_start,offset_end
56      nsph = offset_start-offset_end
57*
58      if (position .gt. 0) then
59        if (lds .lt. ntr) call errquit (
60     &      ' lds < ntr in ecp_cstrans',99, BASIS_ERR)
61        if (ldc .lt. ntr) call errquit (
62     &      ' ldc < ntr in ecp_cstrans',99, BASIS_ERR)
63      else
64        if (lds .lt. nsph) call errquit (
65     &      ' lds < nsph in ecp_cstrans',99, BASIS_ERR)
66        if (ldc .lt. ncart) call errquit (
67     &      ' ldc < ncart in ecp_cstrans',99, BASIS_ERR)
68      end if
69*
70      i = lstart(0,l)+ncart*(ncart-offset_start)
71      j = lstart(1,l)+ncart-offset_start
72C      write (6,*) lstart(0,l),lstart(1,l)
73C      write (6,*) 'l,i,j',l,i,j
74C      call ecp_matpr(csco(i),1,ncart,1,ncart,1,ncart,1,nsph,
75C     &    'Cart to sph','F',81,5)
76C      call ecp_matpr(csco(j),1,ncart,1,ncart,1,nsph,1,ncart,
77C     &    'Sph to cart','F',81,5)
78*
79*     Cartesian to spherical
80*
81      if (indicator .gt. 0) then
82        if (position .gt. 0) then
83          if (mode .gt. 0) then
84            call dgemm ('N','N',ntr,nsph,ncart,one,cart,ldc,
85     &          csco(i),ncart,zero,sph,lds)
86          else
87            call dgemm ('N','T',ntr,nsph,ncart,one,cart,ntr,
88     &          csco(j),ncart,zero,sph,lds)
89          end if
90        else
91          if (mode .gt. 0) then
92            call dgemm ('T','N',nsph,ntr,ncart,one,csco(i),
93     &          ncart,cart,ldc,zero,sph,lds)
94          else
95            call dgemm ('N','N',nsph,ntr,ncart,one,csco(j),
96     &          ncart,cart,ldc,zero,sph,lds)
97          end if
98        end if
99*
100*     Spherical to cartesian
101*
102      else
103        if (position .gt. 0) then
104          if (mode .gt. 0) then
105            call dgemm ('N','N',ntr,ncart,nsph,one,sph,lds,
106     &          csco(j),ncart,zero,cart,ldc)
107          else
108            call dgemm ('N','T',ntr,ncart,nsph,one,sph,lds,
109     &          csco(i),ncart,zero,cart,ldc)
110          end if
111        else
112          if (mode .gt. 0) then
113            call dgemm ('T','N',ncart,ntr,nsph,one,csco(j),
114     &          ncart,sph,lds,zero,cart,ldc)
115          else
116            call dgemm ('N','N',ncart,ntr,nsph,one,csco(i),
117     &          ncart,sph,lds,zero,cart,ldc)
118          end if
119        end if
120      end if
121*
122      return
123      end
124