1C> \ingroup wfn1 2C> @{ 3C> 4C> \brief Evaluate the power of a matrix 5C> 6C> This subroutine calculates \f$A' = A^{p}\f$ where \f$A\f$ is a 7C> square matrix and \f$p\f$ is a floating point number. The function 8C> evaluates the expression in place. I.e. the input matrix is 9C> overwritten with the result. 10C> 11 subroutine wfn1_mat_power(n,a,lda,s,lds,p) 12 implicit none 13c 14#include "mafdecls.fh" 15#include "errquit.fh" 16c 17 integer n !< [Input] The dimension of the matrix 18 integer lda !< [Input] The leading dimension of array \f$A\f$. 19 integer lds !< [Input] The leading dimension of array \f$S\f$. 20c 21 double precision a(lda,n) !< [In/Output] The matrix 22 !< - On input the matrix to calculate the power off 23 !< - On output the input matrix to the power \f$p\f$ 24 double precision s(lds,n) !< [Input] The metric matrix 25 double precision p !< [Input] The power 26c 27 integer l_matvec, k_matvec 28 integer l_diag, k_diag 29 integer l_metric, k_metric 30 integer l_work, k_work 31 integer lwork, info 32 integer itype 33 integer i, j, k 34 double precision size 35c 36 return ! LAPACK dsygv does not work! 37c 38 itype = 1 39 if (.not.ma_push_get(MT_DBL,n*n,"matvec",l_matvec,k_matvec)) 40 & call errquit("wfn1_mat_power: could not allocate matvec",0, 41 & MA_ERR) 42 if (.not.ma_push_get(MT_DBL,n*n,"metric",l_metric,k_metric)) 43 & call errquit("wfn1_mat_power: could not allocate metric",0, 44 & MA_ERR) 45 if (.not.ma_push_get(MT_DBL,n,"diag",l_diag,k_diag)) 46 & call errquit("wfn1_mat_power: could not allocate diag",0, 47 & MA_ERR) 48c 49 do j = 1, n 50 do i = 1, n 51 dbl_mb(k_matvec+(j-1)*n+i-1) = 0.0d0 52 dbl_mb(k_metric+(j-1)*n+i-1) = 0.0d0 53 enddo 54 dbl_mb(k_diag+j-1) = 0.0d0 55 enddo 56 do j = 1, n 57 do i = 1, j 58 dbl_mb(k_matvec+(j-1)*n+i-1) = a(i,j) 59 dbl_mb(k_metric+(j-1)*n+i-1) = s(i,j) 60 enddo 61 enddo 62c 63 info = 0 64 lwork = -1 65 call dsygv(itype,'V','U',n,dbl_mb(k_matvec),n,dbl_mb(k_metric),n, 66 & dbl_mb(k_diag),size,lwork,info) 67 if (info.ne.0) call errquit("wfn1_mat_power w: info non-zero", 68 & info,UERR) 69 lwork = size+0.5d0 70cDEBUG 71c lwork = 3*n-1 72 write(*,*)'wfn1_mat_power: lwork, n: ',lwork,n,size 73cDEBUG 74 if (.not.ma_push_get(MT_DBL,lwork,"work",l_work,k_work)) 75 & call errquit("wfn1_mat_power: could not allocate work",0, 76 & MA_ERR) 77 info = 0 78 itype = 1 79 call dsygv(itype,'V','U',n,dbl_mb(k_matvec),n,dbl_mb(k_metric),n, 80 & dbl_mb(k_diag),dbl_mb(k_work),lwork,info) 81 if (info.ne.0) call errquit("wfn1_mat_power: info non-zero", 82 & info,UERR) 83cDEBUG 84 do i = 0, n-1 85 write(*,*)'wfn1_mat_power: diag: ',i+1,dbl_mb(k_diag+i) 86 enddo 87 write(*,*)'wfn1_mat_power: vectors' 88 call hess_hssout(dbl_mb(k_matvec),n,n,n) 89cDEBUG 90 do i = 0, n-1 91c dbl_mb(k_diag+i) = max(dbl_mb(k_diag+i),0.0d0)**p 92 dbl_mb(k_diag+i) = max(dbl_mb(k_diag+i),0.0d0) 93 enddo 94 do j = 1, n 95 do i = 1, n 96 a(i,j) = 0.0d0 97 do k = 1, n 98 a(i,j) = a(i,j) + dbl_mb(k_diag+k-1) 99 & * dbl_mb(k_matvec+(k-1)*n+i-1) 100 & * dbl_mb(k_matvec+(k-1)*n+j-1) 101 enddo ! k 102 enddo ! i 103 enddo ! j 104 if (.not.ma_pop_stack(l_work)) 105 & call errquit("wfn1_mat_power: could not deallocate work", 106 & 0,MA_ERR) 107 if (.not.ma_pop_stack(l_diag)) 108 & call errquit("wfn1_mat_power: could not deallocate diag", 109 & 0,MA_ERR) 110 if (.not.ma_pop_stack(l_metric)) 111 & call errquit("wfn1_mat_power: could not deallocate metric", 112 & 0,MA_ERR) 113 if (.not.ma_pop_stack(l_matvec)) 114 & call errquit("wfn1_mat_power: could not deallocate matvec", 115 & 0,MA_ERR) 116 end 117C> 118C> @} 119