Line data    Source code

       1              : 
       2              : #if defined HAVE_CONFIG_H
       3              : #include "config.h"
       4              : #endif
       5              : 
       6              : #include "abi_common.h"
       7              : 
       8              : module m_tdep_constraints
       9              : 
      10              :   use defs_basis
      11              :   use m_errors
      12              :   use m_xmpi
      13              : 
      14              :   implicit none
      15              : 
      16              :   type S_product
      17              :     double precision, allocatable :: SS  (:,:,:)
      18              :     double precision, allocatable :: SSS (:,:,:,:)
      19              :     double precision, allocatable :: SSSS(:,:,:,:,:)
      20              :   end type S_product
      21              : 
      22              :   type Asr_Rot
      23              :     double precision, allocatable :: ABG (:,:,:)
      24              :     double precision, allocatable :: ABGD(:,:,:,:)
      25              :     double precision, allocatable :: ABGDE(:,:,:,:,:)
      26              :   end type Asr_Rot
      27              : 
      28              :   type,public :: Constraints_type
      29              :     type(S_product),allocatable :: Sprod(:,:)
      30              :     type(Asr_Rot),allocatable :: AsrRot3(:,:,:)
      31              :     type(Asr_Rot),allocatable :: AsrRot4(:,:,:,:)
      32              :   end type Constraints_type
      33              : 
      34              :   public :: tdep_calc_orthonorm
      35              : 
      36              : contains
      37              : 
      38              : !====================================================================================================
      39              : 
      40          204 : subroutine tdep_calc_orthonorm(dim1,dim2,nindep,vect)
      41              : 
      42              :   integer, intent(in) :: dim1,dim2
      43              :   integer, intent(out) :: nindep
      44              :   double precision, intent(inout) :: vect(dim1,dim2)
      45              : 
      46              :   integer :: ii,jj,kk
      47              :   double precision :: prod_scal
      48              : 
      49              : ! Filter non-zero vectors
      50          204 :   ii=0
      51       539556 :   do kk=1,dim2
      52     51280710 :     if (sum(abs(vect(:,kk))).gt.tol8) then
      53        10106 :       ii=ii+1
      54      1462868 :       vect(:,ii)=vect(:,kk)
      55              :     end if
      56              :   end do
      57          204 :   nindep=ii
      58              : 
      59              : ! Gram-Schmidt orthogonalization
      60        10252 :   do kk=2,nindep
      61      3502333 :     do jj=1,kk-1
      62    478515426 :       prod_scal=sum(vect(:,jj)*vect(:,jj))
      63      3502129 :       if (abs(prod_scal).gt.tol8) then
      64     40441616 :         vect(:,kk)=vect(:,kk)-sum(vect(:,kk)*vect(:,jj))/prod_scal*vect(:,jj)
      65              :       end if
      66              :     end do
      67              :   end do
      68              : 
      69              : ! Store the non-zero vectors and normalize
      70              :   ii=0
      71        10310 :   do kk=1,nindep
      72      1462868 :     prod_scal=sum(vect(:,kk)*vect(:,kk))
      73        10310 :     if (abs(prod_scal).gt.tol8) then
      74          370 :       ii=ii+1
      75        88695 :       vect(:,ii)=vect(:,kk)/dsqrt(prod_scal)
      76              :     end if
      77              :   end do
      78       529450 :   do kk=nindep+1,dim2
      79     49817842 :     vect(:,kk)=zero
      80              :   end do
      81          204 :   nindep=ii
      82              : 
      83          204 : end subroutine tdep_calc_orthonorm
      84              : 
      85              : !====================================================================================================
      86              : 
      87            0 : end module m_tdep_constraints