1*
2* $Id$
3*
4
5*     ***********************************
6*     *					*
7*     *	       c_geodesic2_init		*
8*     *					*
9*     ***********************************
10*
11*     Uses - geodesic common block
12*
13
14      subroutine c_geodesic2_init()
15      implicit none
16
17#include "bafdecls.fh"
18#include "c_geodesic2.fh"
19#include "errquit.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*     *		c_geodesic2_finalize	*
58*     *					*
59*     ***********************************
60*
61*     Uses - c_geodesic2 common block
62*
63
64      subroutine g2eodesic2_finalize()
65      implicit none
66
67#include "bafdecls.fh"
68#include "c_geodesic2.fh"
69#include "errquit.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('c_geodesic2_finalize:freeing heap memory',0,
86     &       MA_ERR)
87
88      return
89      end
90
91
92
93*     ***********************************
94*     *					*
95*     *	      c_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 - c_geodesic2 common block
105*
106      subroutine c_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 "c_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
146*     **** calculate A=<Y|H> ****
147      call Dneall_ffm_Multiply(0,Y,H,npack1,dbl_mb(A(1)))
148
149*     **** calculate Q=(I-YYt)H - should not be necessary but just in case ****
150      call Dneall_fmf_Multiply(0,Y,npack1,
151     >                         dbl_mb(A(1)),1.0d0,
152     >                         dcpl_mb(Q(1)),0.0d0)
153      call daxpy(2*npack1*nemax,(-1.0d0),H,1,dcpl_mb(Q(1)),1)
154      call dscal(2*npack1*nemax,(-1.0d0),dcpl_mb(Q(1)),1)
155
156*     **** calculate QR using Modified Gram-Schmidt ****
157      call Dneall_fm_QR(0,dcpl_mb(Q(1)),npack1,dbl_mb(R(1)))
158
159*     **** generate T from A and R ****
160*       -     -
161*  T = |A, -R^t|
162*      |R,  0  |
163*       -     -
164      call Dneall_AR_to_4m(0,dbl_mb(A(1)),dbl_mb(R(1)),dbl_mb(T(1)))
165
166*     **** Factor T--> V,W,and S ****
167      call Dneall_4m_FactorSkew(0, dbl_mb(T(1)),
168     >                             dbl_mb(V(1)),
169     >                             dbl_mb(W(1)),
170     >                             dbl_mb(S(1)))
171
172*     **** calculate dE ****
173      dE = 2.0d0*c_electron_eorbit_noocc(H)
174
175
176*     **** deallocate tmp space ****
177      value = Dneall_m_free(T)
178      if (.not. value)
179     > call errquit('geodesic2_start:popping stack',1, MA_ERR)
180
181      call nwpw_timing_end(10)
182      return
183      end
184
185*     ***********************************
186*     *					*
187*     *	      c_geodesic2_generate_T	*
188*     *					*
189*     ***********************************
190*
191*     This routine determines T.  T is defined
192* to be a 2nx2n skew symmetric matrix.
193*
194*       -     -
195*  T = |A, -R^t|
196*      |R,  0  |
197*       -     -
198*
199*     Entry - n: dimension of matrices A and R
200*             A: an nxn skew symmetric matrix
201*             R: an nxn matrix
202*     Exit - T: a 2nx2n skew symmetric matrix
203
204      subroutine geodesic2_generate_T(n,A,R,T)
205      implicit none
206      integer n
207      real*8 A(n,n)
208      real*8 R(n,n)
209      real*8 T(2*n,2*n)
210
211*     **** local variables ****
212      integer i,j
213
214      call dcopy(4*n*n,0.0d0,0,T,1)
215
216*     **** copy A to upper-left of T ****
217      do j=1,n
218      do i=1,n
219         T(i,j) = A(i,j)
220      end do
221      end do
222
223*     **** copy R to lower-left of T ****
224      do j=1,n
225      do i=1,n
226         T(i+n,j) = R(i,j)
227      end do
228      end do
229
230*     **** copy -R^t to upper-right of T ****
231      do j=1,n
232      do i=1,n
233         T(i,j+n) = -R(j,i)
234      end do
235      end do
236
237      return
238      end
239
240
241*     ***********************************
242*     *					*
243*     *		c_geodesic2_get		*
244*     *					*
245*     ***********************************
246*
247*     This routine calculates
248*
249*   Ynew = Yold*M(t) + Q*N(t)
250*
251*   where
252*        -    -               - -
253*       | M(t) | = Exp(t*T)* | I |
254*       | N(t) |             | 0 |
255*        -    -               - -
256*
257      subroutine c_geodesic2_get(t,Yold,Ynew)
258      implicit none
259      real*8     t
260      complex*16 Yold(*)
261      complex*16 Ynew(*)
262
263#include "bafdecls.fh"
264#include "errquit.fh"
265#include "geodesic2.fh"
266
267*     **** local variables ****
268      logical value
269      integer npack1,ispin,ne(2)
270      integer ms,n,j,k,shift1,shift2
271      integer MM(2),NN(2)
272
273*     **** external functions ****
274      logical  Dneall_m_push_get,Dneall_m_pop_stack
275      external Dneall_m_push_get,Dneall_m_pop_stack
276
277
278      call Pack_npack(1,npack1)
279
280*     **** allocate tmp space ****
281      value =           Dneall_m_push_get(0,MM)
282      value = value.and.Dneall_m_push_get(0,NN)
283      if (.not. value)
284     >   call errquit('geodesic2_get: pushing stack',0, MA_ERR)
285
286      call geodesic2_get_MandN(t,dbl_mb(MM(1)),dbl_mb(NN(1)))
287
288      call Dneall_fmf_Multiply(0,
289     >                         Yold,npack1,
290     >                         dbl_mb(MM(1)),
291     >                         1.0d0,
292     >                         Ynew,
293     >                         0.0d0)
294      call Dneall_fmf_Multiply(0,
295     >                         dcpl_mb(Q(1)),npack1,
296     >                         dbl_mb(NN(1)),
297     >                         1.0d0,
298     >                         Ynew,
299     >                         1.0d0)
300
301*     **** deallocate tmp space ****
302      value =           Dneall_m_pop_stack(NN)
303      value = value.and.Dneall_m_pop_stack(MM)
304      if (.not. value)
305     > call errquit('geodesic2_get:popping stack',1, MA_ERR)
306
307      return
308      end
309
310
311*     ***********************************
312*     *					*
313*     *		c_geodesic2_transport	*
314*     *					*
315*     ***********************************
316*
317*     This routine calculates
318*
319*   Hnew = Hold*M(t)    + Yold*R^t*N(t)
320*
321*   where
322*        -    -               - -
323*       | M(t) | = Exp(t*T)* | I |
324*       | N(t) |             | 0 |
325*        -    -               - -
326*
327      subroutine c_geodesic2_transport(t,Yold,Hnew)
328      implicit none
329      real*8     t
330      complex*16 Yold(*)
331      complex*16 Hnew(*)
332
333#include "bafdecls.fh"
334#include "c_geodesic2.fh"
335#include "errquit.fh"
336
337
338*     **** local variables ****
339      logical value
340      integer npack1
341      integer MM(2),NN(2),TT(2)
342
343*     **** external functions ****
344      logical  Dneall_m_push_get,Dneall_m_pop_stack
345      external Dneall_m_push_get,Dneall_m_pop_stack
346
347
348      call Pack_npack(1,npack1)
349
350*     **** allocate tmp space ****
351      value =           Dneall_m_push_get(0,MM)
352      value = value.and.Dneall_m_push_get(0,NN)
353      value = value.and.Dneall_m_push_get(0,TT)
354      if (.not. value)
355     >   call errquit('geodesic2_transport: pushing stack',0, MA_ERR)
356
357      call geodesic2_get_MandN(t,dbl_mb(MM(1)),dbl_mb(NN(1)))
358
359*     **** TT(t) = -R^t*NN(t) ****
360      call Dneall_mmm_Multiply2ab(0,
361     >                          dbl_mb(R(1)),
362     >                          dbl_mb(NN(1)),
363     >                          (-1.0d0),
364     >                          dbl_mb(TT(1)),
365     >                          0.0d0)
366
367
368*     *** Hnew <-- Hold*M(t) + Yold*TT(t) ***
369      call Dneall_fmf_Multiply(0,
370     >                         dcpl_mb(Hold(1)),npack1,
371     >                         dbl_mb(MM(1)),
372     >                         1.0d0,
373     >                         Hnew,
374     >                         0.0d0)
375      call Dneall_fmf_Multiply(0,
376     >                         Yold,npack1,
377     >                         dbl_mb(TT(1)),
378     >                         1.0d0,
379     >                         Hnew,
380     >                         1.0d0)
381
382*     **** deallocate tmp space ****
383      value =           Dneall_m_pop_stack(TT)
384      value = value.and.Dneall_m_pop_stack(NN)
385      value = value.and.Dneall_m_pop_stack(MM)
386      if (.not. value)
387     > call errquit('c_geodesic2_transport:popping stack',1,MA_ERR)
388
389
390      return
391      end
392
393      subroutine c_geodesic2_checkMN(name,n,MM,NN)
394      implicit none
395      character*(*) name
396      integer n
397      real*8 MM(n,n),NN(n,n)
398
399*     **** local variables ****
400      integer i,j,k
401      real*8 tmp(n,n),sum
402
403      do j=1,n
404      do i=1,n
405        sum = 0.0d0
406        do k=1,n
407          sum = sum + MM(k,i)*NN(k,j)
408        end do
409        tmp(i,j) = sum
410      end do
411      end do
412
413      write(*,*)
414      write(*,*) "checkMN:",name
415      do i=1,n
416        write(*,'(22F8.3)') (tmp(i,j),j=1,n)
417      end do
418      write(*,*)
419      return
420      end
421
422*     ***********************************
423*     *					*
424*     *	      c_geodesic2_get_MandN	*
425*     *					*
426*     ***********************************
427*
428*   This routine returns
429*        -    -               - -
430*       | M(t) | = Exp(t*T)* | I |
431*       | N(t) |             | 0 |
432*        -    -               - -
433*   where
434*
435*      T =  U*Sigma*U^H, with U=(V+iW)
436*
437*      is a skew matrix that is decomposed into V,W,and Sigma
438*
439      subroutine c_geodesic2_get_MandN(t,M,N)
440      implicit none
441      real*8 t
442      real*8 M(*),N(*)
443
444#include "bafdecls.fh"
445#include "c_geodesic2.fh"
446#include "errquit.fh"
447
448
449*     **** local variables ****
450      logical value
451      integer ispin,ne(2)
452      integer ms,i,j,nn,shift,shift2,shift3
453      integer RR(2)
454
455*     **** external functions ****
456      logical  Dneall_4m_push_get,Dneall_m_pop_stack
457      external Dneall_4m_push_get,Dneall_m_pop_stack
458
459*     **** allocate tmp space ****
460      value = Dneall_4m_push_get(0,RR)
461      if (.not. value)
462     >   call errquit('c_geodesic2_get_MandN: pushing stack',0, MA_ERR)
463
464      call Dneall_4m_RotateSkew(0,t,
465     >                          dbl_mb(V(1)),
466     >                          dbl_mb(W(1)),
467     >                          dbl_mb(S(1)),
468     >                          dbl_mb(RR(1)) )
469
470      call Dneall_4m_to_MN(0,dbl_mb(RR(1)),M,N)
471
472*     **** deallocate tmp space ****
473      value = Dneall_m_pop_stack(RR)
474      if (.not. value)
475     > call errquit('c_geodesic2_get_MandN:popping stack',1, MA_ERR)
476
477      return
478      end
479