1*
2* $Id$
3*
4
5*     ***********************************
6*     *					*
7*     *		geodesic2_init		*
8*     *					*
9*     ***********************************
10*
11*     Uses - geodesic common block
12*
13
14      subroutine geodesic2_init()
15      implicit none
16#include "errquit.fh"
17
18#include "bafdecls.fh"
19#include "geodesic2.fh"
20
21*     **** local variables ****
22      logical value
23      integer npack1
24
25*     **** external functions ****
26      integer  psi_ne,psi_neq
27      external psi_ne,psi_neq
28      logical  Dneall_m_allocate,Dneall_4m_allocate
29      external Dneall_m_allocate,Dneall_4m_allocate
30
31      call Pack_npack(1,npack1)
32c      nemax = psi_neq(1)+psi_neq(2)
33c      ne(1) = psi_neq(1)
34c      ne(2) = psi_neq(2)
35
36      value = BA_alloc_get(mt_dcpl,npack1*(psi_neq(1)+psi_neq(2)),
37     >                     'Hold',Hold(2),Hold(1))
38      value = value.and.
39     >        BA_alloc_get(mt_dcpl,npack1*(psi_neq(1)+psi_neq(2)),
40     >                     'Q',Q(2),Q(1))
41      value = value.and.Dneall_m_allocate(0,R)
42      value = value.and.Dneall_m_allocate(0,A)
43      value = value.and.Dneall_4m_allocate(0,V)
44      value = value.and.Dneall_4m_allocate(0,W)
45
46      value = value.and.
47     >        BA_alloc_get(mt_dbl,2*(psi_ne(1)+psi_ne(2)),
48     >                     'S',S(2),S(1))
49      if (.not. value) call errquit('geodesic2_init:allocating heap',0,
50     &       MA_ERR)
51
52      return
53      end
54
55*     ***********************************
56*     *					*
57*     *		geodesic2_finalize	*
58*     *					*
59*     ***********************************
60*
61*     Uses - geodesic common block
62*
63
64      subroutine geodesic2_finalize()
65      implicit none
66#include "errquit.fh"
67
68#include "bafdecls.fh"
69#include "geodesic2.fh"
70
71*     **** local variables ****
72      logical value
73
74      logical  Dneall_m_free
75      external Dneall_m_free
76
77      value =           BA_free_heap(S(2))
78      value = value.and.Dneall_m_free(W)
79      value = value.and.Dneall_m_free(V)
80      value = value.and.Dneall_m_free(A)
81      value = value.and.Dneall_m_free(R)
82      value = value.and.BA_free_heap(Q(2))
83      value = value.and.BA_free_heap(Hold(2))
84      if (.not. value)
85     >  call errquit('geodesic2_finalize:freeing heap memory',0,
86     &       MA_ERR)
87
88      return
89      end
90
91
92
93*     ***********************************
94*     *					*
95*     *		geodesic2_start		*
96*     *					*
97*     ***********************************
98*
99*     This routine determines the pxp matrices R and YA, and
100* the orthogonal nxp matrix Q.   Q and R are determined from
101* the QR decomposition of the projected direction (I-YY^t)H, and
102* YH is defined as the Lagrange Multiplier pxp matrix Y^tH.
103*
104*     Uses - geodesic2 common block
105*
106      subroutine geodesic2_start(Y,H,max_sigma,dE)
107      implicit none
108      complex*16 Y(*)
109      complex*16 H(*)
110      real*8     max_sigma,dE
111
112#include "bafdecls.fh"
113#include "errquit.fh"
114#include "geodesic2.fh"
115
116*     **** local variables ****
117      logical value
118      integer npack1,nemax
119      integer ms,n,ispin,ne(2)
120      integer shift,shift2,i,j
121      integer T(2)
122
123*     **** external functions ****
124      logical  Dneall_4m_allocate,Dneall_m_free
125      integer  psi_ispin,psi_neq
126      real*8   electron_eorbit_noocc
127      external Dneall_4m_allocate,Dneall_m_free
128      external psi_ispin,psi_neq
129      external electron_eorbit_noocc
130
131      call nwpw_timing_start(10)
132      call Pack_npack(1,npack1)
133      nemax = psi_neq(1) + psi_neq(2)
134      ispin = psi_ispin()
135      ne(1) = psi_neq(1)
136      ne(2) = psi_neq(2)
137
138*     **** allocate tmp space ****
139      value = Dneall_4m_allocate(0,T)
140      if (.not. value)
141     >   call errquit('geodesic2_start: pushing stack',0, MA_ERR)
142
143*     **** Hold <-- H ****
144      !call dcopy(2*npack1*nemax,H,1,dcpl_mb(Hold(1)),1)
145      call Parallel_shared_vector_copy(.true.,2*npack1*nemax,
146     >                                 H,dcpl_mb(Hold(1)))
147
148*     **** calculate A=<Y|H> ****
149      call Dneall_ffm_Multiply(0,Y,H,npack1,dbl_mb(A(1)))
150
151*     **** calculate Q=(I-YYt)H - should not be necessary but just in case ****
152      call Dneall_fmf_Multiply(0,Y,npack1,
153     >                         dbl_mb(A(1)),1.0d0,
154     >                         dcpl_mb(Q(1)),0.0d0)
155      call DAXPY_OMP(2*npack1*nemax,(-1.0d0),H,1,dcpl_mb(Q(1)),1)
156      call DSCAL_OMP(2*npack1*nemax,(-1.0d0),dcpl_mb(Q(1)),1)
157
158*     **** calculate QR using Modified Gram-Schmidt ****
159      call Dneall_fm_QR(0,dcpl_mb(Q(1)),npack1,dbl_mb(R(1)))
160
161*     **** generate T from A and R ****
162*       -     -
163*  T = |A, -R^t|
164*      |R,  0  |
165*       -     -
166      call Dneall_AR_to_4m(0,dbl_mb(A(1)),dbl_mb(R(1)),dbl_mb(T(1)))
167
168*     **** Factor T--> V,W,and S ****
169      call Dneall_4m_FactorSkew(0, dbl_mb(T(1)),
170     >                             dbl_mb(V(1)),
171     >                             dbl_mb(W(1)),
172     >                             dbl_mb(S(1)))
173
174*     **** calculate dE ****
175      dE = 2.0d0*electron_eorbit_noocc(H)
176
177
178*     **** deallocate tmp space ****
179      value = Dneall_m_free(T)
180      if (.not. value)
181     > call errquit('geodesic2_start:popping stack',1, MA_ERR)
182
183      call nwpw_timing_end(10)
184      return
185      end
186
187*     ***********************************
188*     *					*
189*     *	      geodesic2_generate_T	*
190*     *					*
191*     ***********************************
192*
193*     This routine determines T.  T is defined
194* to be a 2nx2n skew symmetric matrix.
195*
196*       -     -
197*  T = |A, -R^t|
198*      |R,  0  |
199*       -     -
200*
201*     Entry - n: dimension of matrices A and R
202*             A: an nxn skew symmetric matrix
203*             R: an nxn matrix
204*     Exit - T: a 2nx2n skew symmetric matrix
205
206      subroutine geodesic2_generate_T(n,A,R,T)
207      implicit none
208      integer n
209      real*8 A(n,n)
210      real*8 R(n,n)
211      real*8 T(2*n,2*n)
212
213*     **** local variables ****
214      integer i,j
215
216      call dcopy(4*n*n,0.0d0,0,T,1)
217
218*     **** copy A to upper-left of T ****
219      do j=1,n
220      do i=1,n
221         T(i,j) = A(i,j)
222      end do
223      end do
224
225*     **** copy R to lower-left of T ****
226      do j=1,n
227      do i=1,n
228         T(i+n,j) = R(i,j)
229      end do
230      end do
231
232*     **** copy -R^t to upper-right of T ****
233      do j=1,n
234      do i=1,n
235         T(i,j+n) = -R(j,i)
236      end do
237      end do
238
239      return
240      end
241
242
243*     ***********************************
244*     *					*
245*     *		geodesic2_get		*
246*     *					*
247*     ***********************************
248*
249*     This routine calculates
250*
251*   Ynew = Yold*M(t) + Q*N(t)
252*
253*   where
254*        -    -               - -
255*       | M(t) | = Exp(t*T)* | I |
256*       | N(t) |             | 0 |
257*        -    -               - -
258*
259      subroutine geodesic2_get(t,Yold,Ynew)
260      implicit none
261      real*8     t
262      complex*16 Yold(*)
263      complex*16 Ynew(*)
264
265#include "bafdecls.fh"
266#include "errquit.fh"
267#include "geodesic2.fh"
268
269*     **** local variables ****
270      logical value
271      integer npack1,ispin,ne(2)
272      integer ms,n,j,k,shift1,shift2
273      integer MM(2),NN(2)
274
275*     **** external functions ****
276      logical  Dneall_m_push_get,Dneall_m_pop_stack
277      external Dneall_m_push_get,Dneall_m_pop_stack
278
279
280      call Pack_npack(1,npack1)
281
282*     **** allocate tmp space ****
283      value =           Dneall_m_push_get(0,MM)
284      value = value.and.Dneall_m_push_get(0,NN)
285      if (.not. value)
286     >   call errquit('geodesic2_get: pushing stack',0, MA_ERR)
287
288      call geodesic2_get_MandN(t,dbl_mb(MM(1)),dbl_mb(NN(1)))
289
290      call Dneall_fmf_Multiply(0,
291     >                         Yold,npack1,
292     >                         dbl_mb(MM(1)),
293     >                         1.0d0,
294     >                         Ynew,
295     >                         0.0d0)
296      call Dneall_fmf_Multiply(0,
297     >                         dcpl_mb(Q(1)),npack1,
298     >                         dbl_mb(NN(1)),
299     >                         1.0d0,
300     >                         Ynew,
301     >                         1.0d0)
302
303*     **** deallocate tmp space ****
304      value =           Dneall_m_pop_stack(NN)
305      value = value.and.Dneall_m_pop_stack(MM)
306      if (.not. value)
307     > call errquit('geodesic2_get:popping stack',1, MA_ERR)
308
309      return
310      end
311
312
313*     ***********************************
314*     *					*
315*     *		geodesic2_transport	*
316*     *					*
317*     ***********************************
318*
319*     This routine calculates
320*
321*   Hnew = Hold*M(t)    + Yold*R^t*N(t)
322*
323*   where
324*        -    -               - -
325*       | M(t) | = Exp(t*T)* | I |
326*       | N(t) |             | 0 |
327*        -    -               - -
328*
329      subroutine geodesic2_transport(t,Yold,Hnew)
330      implicit none
331#include "errquit.fh"
332      real*8     t
333      complex*16 Yold(*)
334      complex*16 Hnew(*)
335
336#include "bafdecls.fh"
337#include "geodesic2.fh"
338
339
340*     **** local variables ****
341      logical value
342      integer npack1
343      integer MM(2),NN(2),TT(2)
344
345*     **** external functions ****
346      logical  Dneall_m_push_get,Dneall_m_pop_stack
347      external Dneall_m_push_get,Dneall_m_pop_stack
348
349
350      call Pack_npack(1,npack1)
351
352*     **** allocate tmp space ****
353      value =           Dneall_m_push_get(0,MM)
354      value = value.and.Dneall_m_push_get(0,NN)
355      value = value.and.Dneall_m_push_get(0,TT)
356      if (.not. value)
357     >   call errquit('geodesic2_transport: pushing stack',0, MA_ERR)
358
359      call geodesic2_get_MandN(t,dbl_mb(MM(1)),dbl_mb(NN(1)))
360
361*     **** TT(t) = -R^t*NN(t) ****
362      call Dneall_mmm_Multiply2ab(0,
363     >                          dbl_mb(R(1)),
364     >                          dbl_mb(NN(1)),
365     >                          (-1.0d0),
366     >                          dbl_mb(TT(1)),
367     >                          0.0d0)
368
369
370*     *** Hnew <-- Hold*M(t) + Yold*TT(t) ***
371      call Dneall_fmf_Multiply(0,
372     >                         dcpl_mb(Hold(1)),npack1,
373     >                         dbl_mb(MM(1)),
374     >                         1.0d0,
375     >                         Hnew,
376     >                         0.0d0)
377      call Dneall_fmf_Multiply(0,
378     >                         Yold,npack1,
379     >                         dbl_mb(TT(1)),
380     >                         1.0d0,
381     >                         Hnew,
382     >                         1.0d0)
383
384*     **** deallocate tmp space ****
385      value =           Dneall_m_pop_stack(TT)
386      value = value.and.Dneall_m_pop_stack(NN)
387      value = value.and.Dneall_m_pop_stack(MM)
388      if (.not. value)
389     > call errquit('geodesic2_transport:popping stack',1, MA_ERR)
390
391
392      return
393      end
394
395      subroutine geodesic2_checkMN(name,n,MM,NN)
396      implicit none
397      character*(*) name
398      integer n
399      real*8 MM(n,n),NN(n,n)
400
401*     **** local variables ****
402      integer i,j,k
403      real*8 tmp(n,n),sum
404
405      do j=1,n
406      do i=1,n
407        sum = 0.0d0
408        do k=1,n
409          sum = sum + MM(k,i)*NN(k,j)
410        end do
411        tmp(i,j) = sum
412      end do
413      end do
414
415      write(*,*)
416      write(*,*) "checkMN:",name
417      do i=1,n
418        write(*,'(22F8.3)') (tmp(i,j),j=1,n)
419      end do
420      write(*,*)
421      return
422      end
423
424*     ***********************************
425*     *					*
426*     *	      geodesic2_get_MandN	*
427*     *					*
428*     ***********************************
429*
430*   This routine returns
431*        -    -               - -
432*       | M(t) | = Exp(t*T)* | I |
433*       | N(t) |             | 0 |
434*        -    -               - -
435*   where
436*
437*      T =  U*Sigma*U^H, with U=(V+iW)
438*
439*      is a skew matrix that is decomposed into V,W,and Sigma
440*
441      subroutine geodesic2_get_MandN(t,M,N)
442      implicit none
443      real*8 t
444      real*8 M(*),N(*)
445
446#include "bafdecls.fh"
447#include "errquit.fh"
448#include "geodesic2.fh"
449
450
451*     **** local variables ****
452      logical value
453      integer ispin,ne(2)
454      integer ms,i,j,nn,shift,shift2,shift3
455      integer RR(2)
456
457*     **** external functions ****
458      logical  Dneall_4m_push_get,Dneall_m_pop_stack
459      external Dneall_4m_push_get,Dneall_m_pop_stack
460
461*     **** allocate tmp space ****
462      value = Dneall_4m_push_get(0,RR)
463      if (.not. value)
464     >   call errquit('geodesic2_get_MandN: pushing stack',0, MA_ERR)
465
466      call Dneall_4m_RotateSkew(0,t,
467     >                          dbl_mb(V(1)),
468     >                          dbl_mb(W(1)),
469     >                          dbl_mb(S(1)),
470     >                          dbl_mb(RR(1)) )
471
472      call Dneall_4m_to_MN(0,dbl_mb(RR(1)),M,N)
473
474*     **** deallocate tmp space ****
475      value = Dneall_m_pop_stack(RR)
476      if (.not. value)
477     > call errquit('geodesic2_get_MandN:popping stack',1, MA_ERR)
478
479      return
480      end
481