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