Remove QR for CSFs because it does not work.

src/utils/linear_algebra.irp.f

@@ -1133,102 +1133,103 @@ subroutine ortho_svd(A,LDA,m,n)
 
 end
 
-subroutine ortho_qr_withB(A,LDA,B,m,n)
+! QR to orthonormalize CSFs does not work :-(
-  implicit none
+!subroutine ortho_qr_withB(A,LDA,B,m,n)
-  BEGIN_DOC
+!  implicit none
-  ! Orthogonalization using Q.R factorization
+!  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)
 !
-  ! A : Overlap Matrix
+!  integer                        :: LWORK, INFO
+!  integer, allocatable           :: jpvt(:)
+!  double precision, allocatable  :: TAU(:), WORK(:)
+!  double precision, allocatable  :: C(:,:)
+!  double precision :: norm
+!  integer :: i,j
 !
-  ! LDA : leftmost dimension of A
+!  allocate (TAU(min(m,n)), WORK(1))
+!  allocate (jpvt(n))
+!  !print *," In function ortho"
+!  B = A
 !
-  ! m : Number of rows of A
+!  jpvt(1:n)=1
 !
-  ! 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 )
+!  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 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))
+!  deallocate(WORK)
-  !call dgemm('N','N',m,n,n,1.0d0,B,LDA,A,LDA,0.0d0,C,LDA)
+!  allocate(WORK(LWORK))
-  !norm = 0.0d0
+!  call dgeqp3(m, n, A, LDA, jpvt, TAU, WORK, LWORK, INFO )
-  !B = 0.0d0
+!  print *,A
-  !!print *,C
+!  print *,jpvt
-  !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
+!  !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")
+!subroutine ortho_qr_csf(A, LDA, B, m, n) bind(C, name="ortho_qr_csf")
-  use iso_c_binding
+!  use iso_c_binding
-  integer(c_int32_t), value      :: LDA
+!  integer(c_int32_t), value      :: LDA
-  integer(c_int32_t), value      :: m
+!  integer(c_int32_t), value      :: m
-  integer(c_int32_t), value      :: n
+!  integer(c_int32_t), value      :: n
-  integer(c_int16_t)             :: A(LDA,n)
+!  integer(c_int16_t)             :: A(LDA,n)
-  integer(c_int16_t)             :: B(LDA,n)
+!  integer(c_int16_t)             :: B(LDA,n)
-  call ortho_qr_withB(A,LDA,B,m,n)
+!  call ortho_qr_withB(A,LDA,B,m,n)
-end subroutine ortho_qr_csf
+!end subroutine ortho_qr_csf
 
 subroutine ortho_qr(A,LDA,m,n)
   implicit none