C HBMDPP    SOURCE    OF166741  26/05/11    21:15:09     12538          

C=======================================================================
* Calcule la derivee du terme DPsiPsi = (Psi*Psi')/(Psi'*Psi)
* ou: Psi = (L(o)In)*X = Psix*X
C=======================================================================

      SUBROUTINE HBMDPP(NT,NDDL,X,PHI,DPsiPsi)

      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)

      REAL*8 X(*),PHI(*),DPsiPsi(NT,*)

      SEGMENT mwrkl
        REAL*8 PsiV, PsiPsi,BB
        REAL*8 LV(NT),PSI(NT)
        REAL*8 PSIX(NT,NT),PEXSI(NT,NT),PEXV(NT,NT),Ax(NT,NT),MUX(NT,NT)
      ENDSEGMENT

C Fonctions BLAS/LAPACK
      REAL*8    DDOT, DNRM2
      EXTERNAL  DDOT, DNRM2

c   Nombre d'harmoniques
      NHBM = (NT/NDDL - 1)/2

c   Initialisation
      SEGINI,mwrkl
      DO I = 1,NT
        DO J = 1,NT
          DPsiPsi(I,J) = 0.D0
c*          PSIX(I,J) = 0.D0
        ENDDO
      ENDDO

c   Calcul des vecteurs PSI,LV
      CALL HBMDVEC(NT,NHBM,NDDL,X,1.D0,PSI)
      CALL HBMDVEC(NT,NHBM,NDDL,PHI,1.D0,LV)

c   Produits scalaires
      PsiV   = DDOT(NT,PSI,1,PHI,1)
      PSiPsi = DDOT(NT,PSI,1,PSI,1)

c   Produits externes
      CALL PREXT(NT,PSI,PSI,PEXSI)
      CALL PREXT(NT,LV,PSI,PEXV)

c   Construction de Psix*psiv
      DO J=2,2*NHBM,2
        BB = 0.5D0*PsiV*J
        j_z1 = NDDL*(1+(J-1))
        j_z2 = NDDL*(1+(J-2))
        DO I=1,NDDL
          PSIX(j_z2+I,j_z1+I) = BB
          PSIX(j_z1+I,j_z2+I) = -BB
        ENDDO
      ENDDO

*   Construction de Ax = I-(2/psipsi)*PEXSI
      r_z = -2.D0 / PsiPsi
      DO J = 1,NT
        DO I = 1,NT
          Ax(I,J) = r_z * PEXSI(I,J)
        ENDDO
      ENDDO
      DO I = 1,NT
        AX(I,I) = AX(I,I) + 1.D0
      ENDDO

c   Produit MUX = Ax*Psix*psiv
      CALL PRMAT(NT,Ax,PSIX,MUX)

c   Construction de DPsiPsi
      r_z = 1.D0 / PsiPsi
      DO J = 1,NT
        DO I = 1,NT
          DPsiPsi(I,J) = (MUX(I,J)-PEXV(I,J))*r_z
        ENDDO
      ENDDO

      SEGSUP,mwrkl

c      return
      END