1C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
2C NAME
3C     GA_MatPow -- Form V^{-1/2}
4C
5C REVISION
6C     $Id$
7C
8C SYNOPSIS
9      Subroutine GA_MatPow(G_V, Pow, MinEval)
10      Implicit NONE
11#include "errquit.fh"
12      Integer G_V
13      Double Precision Pow, MinEval
14C
15C ARGUMENTS
16C     G_V      Global array handle to the matrix of interest [INOUT]
17C     Pow      Power to which matrix is to be taken [IN]
18C     MinEval  Minimum eigenvalue of input matrix to be retained in
19C              the case of inversion [IN]
20C
21C DESCRIPTION
22C     Forms V^{1/2}, V^{-1/2} or V^{-1} depending on the value of Pow.
23C     Uses the eigenvalue decomposition approach.  In the case of
24C     Pow < 0, MinEval allows filtering of small, potentially
25C     problemmatic eigenvalues.
26C
27C     The input GA is overwritten with the exponentiated result.  It is
28C     NOT guaranteed that the same handle will be returned -- if it is
29C     most efficient, the original V may be destroyed and a new GA created
30C     to hold the result.
31C
32C     This routine responds to the print directive "ga_matpow details",
33C     which is at the Print_Debug level by default.  The print control
34C     context is that of the caller.
35C
36C MEMORY USE
37C     Uses a GA the size of V, and a local array the size of the number
38C     of rows of V.  The eigensolver requires additional memory.
39C
40C     Due to the use of a generalized eigensolver, an additional GA the
41C     size of V is also used.
42C
43C INCLUDE FILES
44#include "global.fh"
45#include "mafdecls.fh"
46#include "numerical_constants.fh"
47#include "stdio.fh"
48#include "util.fh"
49C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
50C CONSTANTS
51      Integer MinChunk
52      Parameter (MinChunk = 64)
53C
54C LOCAL VARIABLES
55      Logical Status, Have_Diag
56      Integer H_Eval, Eval_Idx
57      Integer D, DiagLo, DiagHi
58      Integer G_EV, G_Tmp
59      Integer I, J, VType, Rows, Cols, ILo, IHi, JLo, JHi
60#if defined(USE_GEN_DIAG)
61      Integer G_Junk
62#endif
63      Double precision Elem
64C
65C     ***************************
66C     * Check inputs for sanity *
67C     ***************************
68C
69      Call GA_Check_Handle( G_V, 'GA_MatPow: Input matrix')
70C
71      If ( Pow .ne. FP_Half .AND. Pow .ne. -FP_Half
72     $   .AND. Pow .ne. -FP_One) Call ErrQuit(
73     $   'GA_MatPow: Don''t know how to deal with this exponent',
74     $   Int( 100 * Pow ) , GA_ERR)
75C
76      If ( Util_Print('ga_matpow details', Print_Debug) ) then
77         If ( GA_NodeID() .eq. 0 )
78     $        Call Util_Print_Centered(LuOut, 'GA_MatPow Input', 40,
79     $        .TRUE.)
80         Call GA_Sync
81         Call GA_Print( G_V )
82      EndIf
83C
84C     *********************************
85C     * Get basic info about input GA *
86C     *********************************
87C
88      Call GA_Inquire( G_V, VType, Rows, Cols)
89      If ( Rows .ne. Cols )
90     $   Call Errquit('GA_MatPow: Matrix must be square', Rows-Cols,
91     &       GA_ERR)
92C
93C     ************************
94C     * Find eigenstuff of V *
95C     ************************
96C     Create a global array for the eigenvectors...
97C     And a local array for the eigenvalues (that's the way the routine
98C     does it!)
99C
100      Status = .TRUE.
101      Status = Status .and. GA_Create( VType, Rows, Cols,
102     $   'GA_MatPow eigenvectors', -MinChunk, -MinChunk, G_EV)
103      Status = Status .and. MA_Push_Get( VType, Rows,
104     $   'GA_MatPow eigenvalues', H_EVal, EVal_Idx)
105C
106      If ( .NOT. Status )
107     $   Call ErrQuit('GA_MatPow: Unable to allocate temporaries',
108     $   (Rows+1)*Cols, GA_ERR)
109C
110C     Now we can solve the stupid thing...
111C
112#if defined(USE_GEN_DIAG)
113C     Create an array for the metric and turn it into the unit matrix
114C
115      If ( .NOT. GA_Create( VType, Rows, Cols, 'Metric',
116     $   -MinChunk, -MinChunk, G_Junk) )
117     $   Call ErrQuit('GA_MatPow: Unable to allocate metric',
118     $   Rows*Cols, GA_ERR)
119C
120C     This is clearly serial, but hopefully soon we should be able to
121C     go back to the standard eigensolver instead of the generalized one.
122C
123      Call GA_Zero( G_Junk)
124      If ( GA_NodeID() .eq. 0) then
125         Do I = 1, Rows
126            Call GA_Put(G_Junk, I, I, I, I, FP_One, 1)
127         EndDo
128      EndIf
129C
130#if defined(PARALLEL_DIAG)
131#ifdef SCALAPACK
132         call ga_pdsygv(g_v, g_junk, g_ev, dbl_mb(eval_idx),0)
133#else
134      Call GA_Diag( G_V, G_Junk, G_EV, Dbl_MB(Eval_Idx))
135#endif
136#else
137      Call GA_Diag_Seq( G_V, G_Junk, G_EV, Dbl_MB(Eval_Idx))
138#endif
139C
140      If ( .NOT. GA_Destroy( G_Junk) )
141     $   Call ErrQuit('GA_MatPow: Unable to free metric', 0, GA_ERR)
142C
143#else /* defined(USE_GEN_DIAG) */
144#if defined(PARALLEL_DIAG)
145#ifdef SCALAPACK
146      call ga_pdsyev(g_v, g_ev, dbl_mb(eval_idx),0)
147#else
148      Call GA_Diag_Std( G_V, G_EV, Dbl_MB(Eval_Idx))
149#endif
150#else
151      Call GA_Diag_Std_Seq( G_V, G_EV, Dbl_MB(Eval_Idx))
152#endif
153#endif /* defined(USE_GEN_DIAG) */
154C
155C     Sometimes people may just want a look at the eigenvalues.
156C
157      If ( Util_Print('ga_matpow eigenvalues', Print_Debug)
158     $   .OR. Util_Print('ga_matpow details', Print_Debug) ) then
159         If ( GA_NodeID() .eq. 0 ) then
160            Call Util_Print_Centered(LuOut, 'GA_MatPow Eigenvalues', 40,
161     $           .TRUE.)
162            Do I = 1, Rows, 4
163               Write (LuOut, '(1X, I5, 4(1X, E15.6) )') I,
164     $            (Dbl_MB(Eval_Idx + J - 1), J=I,Min(I+3, Rows) )
165            EndDo
166         EndIf
167         Call GA_Sync
168      EndIf
169      If ( Util_Print('ga_matpow details', Print_Debug) ) then
170         If ( GA_NodeID() .eq. 0 ) then
171            Call Util_Print_Centered(LuOut, 'GA_MatPow Vectors', 40,
172     $           .TRUE.)
173         EndIf
174         Call GA_Sync
175         Call GA_Print( G_EV )
176      EndIf
177C
178C     **************************************************
179C     * Filter out potentially problematic eigenvalues *
180C     **************************************************
181C     Everybody find out what panel of V they hold, and what its
182C     diagonal may be.  Then we can have each node work on
183C     what it "needs"
184C
185      Call GA_Distribution( G_V, GA_NodeID(), Ilo, Ihi,
186     $   Jlo, Jhi)
187C
188C     Consider the diagonal of a matrix.  It divides the upper triangle,
189C     where I (row) < J (column) from the lower triangle, where I > J.
190C
191C     If IHi-JLo > 0 the lower-left corner of the patch is below the diagonal.
192C     If ILo-JHi < 0 the upper-right corner of the patch is above the diagonal.
193C     If the diagonal passes through the patch, these two corners must
194C     be on opposite sides of the diagonal, so the two differences must be
195C     of opposite sign, or their product must be negative.
196C
197C     If the diagonal hits precisely on the ll or ur corner, the
198C     corresponding difference will be zero, and the product will
199C     be zero.  This we must also catch.
200C
201      Have_Diag = ( IHi-JLo ) * ( ILo-JHi ) .le. 0
202C
203      DiagLo = Max( ILo, JLo )
204      DiagHi = Min( IHi, JHi )
205C
206      If ( Have_Diag ) then
207         Do D = DiagLo, DiagHi
208            Elem = Dbl_MB( Eval_Idx + D - 1)
209C
210C           It seems that in the RI-MP2 application, very large
211C           fitting basis sets can lead to a fair number of eigenvalues
212C           which are negative, but small (< ~1e-5).  We suppose that
213C           these arise from small numerical inaccuracies in the
214C           integrals, and perhaps in the eigensolver.  In any case,
215C           requiring strict positive definiteness seems too harsh
216C           in practice.
217C
218c$$$C
219c$$$            If ( Elem .lt. FP_Zero ) then
220c$$$               Write (LuErr, *) 'V is not positive definite! ',
221c$$$     $            D, Elem
222c$$$               GA_MatPow = .FALSE.
223c$$$               Return
224c$$$            EndIf
225C
226            If ( Pow .lt. FP_Zero .AND. Elem .lt. MinEval
227     $           .OR. Elem .lt. FP_Zero) then
228               Write (LuErr, *)
229     $            'Dropping small eigenvalue ', D, Elem
230               Elem = FP_Zero
231            ElseIf (Pow .eq. -FP_Half) then
232               Elem = FP_One / Sqrt( Elem )
233            ElseIf (Pow .eq. -FP_One) then
234               Elem = FP_One / Elem
235            ElseIf (Pow .eq. FP_Half) then
236               Elem = Sqrt( Elem )
237            EndIf
238C
239            Dbl_MB( Eval_Idx + D - 1) = Elem
240C
241         EndDo
242      EndIf
243C
244C     Turn V into a diagonal matrix by copying in our bit of the
245C     diagonal.
246C
247      Call GA_Zero( G_V)
248C
249      If ( Have_Diag ) then
250         Do D = DiagLo, DiagHi
251            Call GA_Put( G_V, D, D, D, D,
252     $         Dbl_MB( Eval_Idx + D - 1), 1)
253         EndDo
254      EndIf
255      Call GA_Sync
256C
257      If ( Util_Print('ga_matpow details', Print_Debug) ) then
258         If ( GA_NodeID() .eq. 0 )
259     $        Call Util_Print_Centered(LuOut, 'GA_MatPow Diagonal', 40,
260     $           .TRUE.)
261         Call GA_Sync
262         Call GA_Print( G_V )
263      EndIf
264C
265C     Now we are done with the local eigenvalue array
266C
267      If ( .NOT. MA_Pop_Stack( H_Eval ) )
268     $   Call ErrQuit('GA_MatPow: Unable to free local eigenvalues', 0,
269     &       MA_ERR)
270C
271C     ****************************
272C     * V^{-1/2} = X' v^{-1/2} X *
273C     ****************************
274C     Note that the eigenvalues, v^{-1/2} is actually stored in the
275C     form of a full matrix.  Currently GA offers no way to go in
276C     and scale rows/columns, which would be more efficient.
277C
278      If ( .NOT. GA_Create( VType, Rows, Cols, 'Temporary matrix',
279     $   -MinChunk, -MinChunk, G_Tmp) )
280     $   Call ErrQuit('GA_MatPow: Unable to allocate temporary matrix',
281     $   Rows*Cols, GA_ERR)
282C
283      Call GA_DGEMM( 'N', 'T', Rows, Cols, Cols, FP_One, G_V, G_EV,
284     $   FP_Zero, G_Tmp)
285      Call GA_DGEMM( 'N', 'N', Rows, Cols, Cols, FP_One, G_EV, G_Tmp,
286     $   FP_Zero, G_V)
287C
288C     Get rid of the global arrays that are not persistent
289C
290      Status = .TRUE.
291      Status = Status .and. GA_Destroy( G_Tmp)
292      Status = Status .and. GA_Destroy( G_EV)
293C
294      If ( .NOT. Status )
295     $   Call ErrQuit('GA_MatPow: Unable to free global arrays', 0,
296     &       GA_ERR)
297C
298      If ( Util_Print('ga_matpow details', Print_Debug) ) then
299         If ( GA_NodeID() .eq. 0 )
300     $        Call Util_Print_Centered(LuOut, 'GA_MatPow Output', 40,
301     $           .TRUE.)
302         Call GA_Sync
303         Call GA_Print( G_V )
304      EndIf
305C
306      Return
307      End
308