LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-10-30 10:18:07 +01:00
Code Releases Activity

Working on LAPACK gramSchmidt. Not working with naive dgeqrf.

Browse Source
This commit is contained in:
v1j4y 2022-07-09 09:00:29 +02:00
parent 0394865bf6
commit bda97b39d6
4 changed files with 136 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
src/csf/cfgCI_interface.f90
View File

@ -46,6 +46,24 @@ module cfunctions
real (kind=C_DOUBLE ),intent(out) :: csftodetmatrix(rowsmax,colsmax)
end subroutine getCSFtoDETTransformationMatrix
end interface
interface
subroutine gramSchmidt(A, m, n, B) bind(C, name='gramSchmidt')
import C_INT32_T, C_INT64_T, C_DOUBLE
integer(kind=C_INT32_T),value,intent(in) :: m
integer(kind=C_INT32_T),value,intent(in) :: n
real (kind=C_DOUBLE ),intent(in) :: A(m,n)
real (kind=C_DOUBLE ),intent(out) :: B(m,n)
end subroutine gramSchmidt
end interface
interface
subroutine gramSchmidt_qp(A, m, n, B) bind(C, name='gramSchmidt_qp')
import C_INT32_T, C_INT64_T, C_DOUBLE
integer(kind=C_INT32_T),value,intent(in) :: m
integer(kind=C_INT32_T),value,intent(in) :: n
real (kind=C_DOUBLE ),intent(in) :: A(m,n)
real (kind=C_DOUBLE ),intent(out) :: B(m,n)
end subroutine gramSchmidt_qp
end interface
end module cfunctions
subroutine f_dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC) &

20
src/csf/cfgCI_utils.c
View File

@ -252,6 +252,26 @@ void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOM
buildTreeDriver(bftree, *NSOMO, MS, NBF);
}
void ortho_qr_csf(double *overlapMatrix, int lda, double *orthoMatrix, int rows, int cols);
void gramSchmidt_qp(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
int i,j;
//for(j=0;j<cols;++j){
// for(i=0;i<rows;++i){
// printf(" %3.2f ",overlapMatrix[j*rows + i]);
// }
// printf("\n");
//}
// Call the function ortho_qr from qp
ortho_qr_csf(overlapMatrix, rows, orthoMatrix, rows, cols);
//for(j=0;j<cols;++j){
// for(i=0;i<rows;++i){
// printf(" %3.2f ",orthoMatrix[j*rows + i]);
// }
// printf("\n");
//}
}
void gramSchmidt(double *overlapMatrix, int rows, int cols, double *orthoMatrix){
// vector

1
src/csf/tree_utils.h
View File

@ -47,6 +47,7 @@ void generateAllBFs(int64_t Isomo, int64_t MS, Tree *bftree, int *NBF, int *NSOM
void getSetBits(int64_t n, int *nsetbits);
void getOverlapMatrix(int64_t Isomo, int64_t MS, double **overlapMatrixptr, int *rows, int *cols, int *NSOMOout);
void getOverlapMatrix_withDet(double *bftodetmatrixI, int rowsbftodetI, int colsbftodetI, int64_t Isomo, int64_t MS, double **overlapMatrixI, int *rowsI, int *colsI, int *NSOMO);
void gramSchmidt_qp(double *overlapMatrix, int rows, int cols, double *orthoMatrix);
void gramSchmidt(double *overlapMatrix, int rows, int cols, double *orthoMatrix);

97
src/utils/linear_algebra.irp.f
View File

@ -1133,6 +1133,103 @@ subroutine ortho_svd(A,LDA,m,n)
end
subroutine ortho_qr_withB(A,LDA,B,m,n)
implicit none
BEGIN_DOC
! Orthogonalization using Q.R factorization
!
! A : Overlap Matrix
!
! LDA : leftmost dimension of A
!
! m : Number of rows of A
!
! n : Number of columns of A
!
! B : Output orthogonal basis
!
END_DOC
integer, intent(in) :: m,n, LDA
double precision, intent(inout) :: A(LDA,n)
double precision, intent(inout) :: B(LDA,n)
integer :: LWORK, INFO
integer, allocatable :: jpvt(:)
double precision, allocatable :: TAU(:), WORK(:)
double precision, allocatable :: C(:,:)
double precision :: norm
integer :: i,j
allocate (TAU(min(m,n)), WORK(1))
allocate (jpvt(n))
!print *," In function ortho"
B = A
jpvt(1:n)=1
LWORK=-1
call dgeqp3( m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
LWORK=max(n,int(WORK(1)))
deallocate(WORK)
allocate(WORK(LWORK))
call dgeqp3(m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
print *,A
print *,jpvt
deallocate(WORK,TAU)
!stop
!LWORK=-1
!call dgeqrf( m, n, A, LDA, TAU, WORK, LWORK, INFO )
!! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
!LWORK=max(n,int(WORK(1)))
!deallocate(WORK)
!allocate(WORK(LWORK))
!call dgeqrf(m, n, A, LDA, TAU, WORK, LWORK, INFO )
!LWORK=-1
!call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
!! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648
!LWORK=max(n,int(WORK(1)))
!deallocate(WORK)
!allocate(WORK(LWORK))
!call dorgqr(m, n, n, A, LDA, TAU, WORK, LWORK, INFO)
!
!allocate(C(LDA,n))
!call dgemm('N','N',m,n,n,1.0d0,B,LDA,A,LDA,0.0d0,C,LDA)
!norm = 0.0d0
!B = 0.0d0
!!print *,C
!do i=1,m
! norm = 0.0d0
! do j=1,n
! norm = norm + C(j,i)*C(j,i)
! end do
! norm = 1.0d0/dsqrt(norm)
! do j=1,n
! B(j,i) = C(j,i)
! end do
!end do
!print *,B
!deallocate(WORK,TAU)
end
subroutine ortho_qr_csf(A, LDA, B, m, n) bind(C, name="ortho_qr_csf")
use iso_c_binding
integer(c_int32_t), value :: LDA
integer(c_int32_t), value :: m
integer(c_int32_t), value :: n
integer(c_int16_t) :: A(LDA,n)
integer(c_int16_t) :: B(LDA,n)
call ortho_qr_withB(A,LDA,B,m,n)
end subroutine ortho_qr_csf
subroutine ortho_qr(A,LDA,m,n)
implicit none
BEGIN_DOC