nl_solv.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
C  |---- version non //------------
!||====================================================================
!||    nl_solv             ../engine/source/implicit/nl_solv.F
!||--- called by ------------------------------------------------------
!||    imp_solv            ../engine/source/implicit/imp_solv.F
!||--- calls      -----------------------------------------------------
!||    al_constraint1_hp   ../engine/source/implicit/nl_solv.F
!||    al_constraint2_hp   ../engine/source/implicit/nl_solv.F
!||    bfgs_0              ../engine/source/implicit/imp_bfgs.F
!||    bfgs_ls             ../engine/source/implicit/imp_bfgs.F
!||    cp_real             ../engine/source/implicit/produt_v.F
!||    crit_ite            ../engine/source/implicit/nl_solv.F
!||    frac_d_hp           ../engine/source/implicit/integrator.F
!||    frac_dd_hp          ../engine/source/implicit/integrator.F
!||    get_max             ../engine/source/implicit/nl_solv.F
!||    imp_stop            ../engine/source/implicit/imp_solv.F
!||    integrator2_hp      ../engine/source/implicit/integrator.F
!||    lin_solv            ../engine/source/implicit/lin_solv.F
!||    line_s              ../engine/source/implicit/nl_solv.F
!||    line_s1             ../engine/source/implicit/nl_solv.F
!||    produt_hp           ../engine/source/implicit/produt_v.F
!||    produt_uhp          ../engine/source/implicit/produt_v.F
!||    produt_uhp2         ../engine/source/implicit/produt_v.F
!||    produt_vmhp         ../engine/source/implicit/produt_v.F
!||    vaxpy_hp            ../engine/source/implicit/produt_v.F
!||--- uses       -----------------------------------------------------
!||    dsgraph_mod         ../engine/share/modules/dsgraph_mod.F
!||    groupdef_mod        ../common_source/modules/groupdef_mod.F
!||    intbufdef_mod       ../common_source/modules/interfaces/intbufdef_mod.F90
!||====================================================================
      SUBROUTINE nl_solv(NDDL  ,IDDL  ,NDOF  ,IKC   ,D     ,
     1                   DR    ,NNZ   ,IADK  ,JDIK  ,DIAG_K,
     2                   LT_K  ,F     ,NDDLI ,IADI  ,JDII  ,   
     3                   DIAG_I,LT_I  ,ITOK  ,IADM  ,JDIM  ,   
     4                   DIAG_M,LT_M  ,R02   ,DD    ,DDR   ,   
     5                   ITASK0,IT    ,ITC   ,RU0   ,ROLD  ,
     6                   IDIV ,INPRINT,ICPREC,ISTOP ,E02   ,
     7                   DE0  ,EIMP   ,INLOC ,NDDL0 ,LS    ,
     8                   U02  ,GAP    ,ITAB  ,FR_ELEM,IAD_ELEM,
     9                   W_DDL,A      ,AR    ,V     ,MS    ,
     A                   X    ,IPARI  ,INTBUF_TAB   ,NUM_IMP,
     B                   NS_IMP,NE_IMP,NSREM ,NSL   ,ICONT ,
     C                   GRAPHE,FAC_K, IPIV_K, NK   ,NMONV ,
     D                   IMONV ,MONVOL,IGRSURF,FR_MV ,
     E                   VOLMON,IBFV  ,SKEW  ,XFRAME,MUMPS_PAR,
     F                   CDDLP ,IND_IMP,NBINTC,INTLIST,NEWFRONT,
     G                   ISENDTO,IRECVFROM,IRBE3 ,LRBE3,NDIV ,
     H                   ICONT0 ,ISIGN ,FEXT  ,DG    ,DGR   ,
     I                   DG0   ,DGR0  ,RFEXT ,LS1   ,NODFT ,
     J                   NODLT ,IRBE2 ,LRBE2 ,IDIV0,
     K                   RELRES,anew_stif )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE dsgraph_mod
      USE intbufdef_mod
      USE groupdef_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#if defined(MUMPS5)
#include "dmumps_struc.h"
#endif
#include "com01_c.inc"
#include "com04_c.inc"
#include "com08_c.inc"
#include "impl1_c.inc"
#include "impl2_c.inc"
#include "units_c.inc"
#include "task_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
C----------resol [K]{X}={F}------ 
      INTEGER  NDDL  ,NNZ ,IADK(*),JDIK(*),IADM(*),JDIM(*),ITASK0,
     .         NDDLI  ,IADI(*),JDII(*),ITOK(*),ICPREC,INLOC(*),
     .         NDOF(*),IDDL(*),IKC(*),IDIV  ,INPRINT,ISTOP,
     .         IT,ITC,NDDL0,ITAB(*),FR_ELEM(*),IAD_ELEM(*),W_DDL(*)
      INTEGER  NE_IMP(*),NSREM ,NSL,ICONT,ISIGN,IMP1,
     .         IPARI(*) ,NUM_IMP(*),NS_IMP(*), IPIV_K(*), NK 
      INTEGER NMONV,IMONV(*),MONVOL(*),FR_MV(*),
     .        IBFV(*),IND_IMP(*),IRBE3(*) ,LRBE3(*),NDIV ,ICONT0,
     .        IRBE2(*),LRBE2(*) 
      INTEGER NEWFRONT(*),NBINTC,INTLIST(*),ISENDTO(*),IRECVFROM(*),
     .        NODFT ,NODLT,IDIV0
C     REAL
      my_real
     .  DIAG_K(*),LT_K(*),F(*),R02,DIAG_M(*),LT_M(*),
     .  DIAG_I(*),LT_I(*),D(3,*),DR(3,*),DD(3,*),DDR(3,*),
     .  RU0,ROLD,E02 ,EIMP,LSTOL,DE0,LS(*),U02,GAP,RBID
      my_real
     .  A(3,*),AR(3,*),V(3,*),X(3,*),MS(*),FAC_K(*),
     .  VOLMON(*),SKEW(*),XFRAME(*),FEXT(*),DG(3,*),DGR(3,*),
     .  DG0(3,*),DGR0(3,*),LS1(*)   ,RFEXT, RELRES(*)
C-----LS(1):R0,LS(2):S0,LS(3):S1 ; LS1(1): EN, LS1(2):SN;LS1(3):DE
      TYPE(prgraph) :: GRAPHE(*)
      INTEGER CDDLP(*)
      TYPE(INTBUF_STRUCT_) INTBUF_TAB(*)
      TYPE (SURF_)   , DIMENSION(NSURF)   :: IGRSURF
C
#ifdef MUMPS5
      TYPE(DMUMPS_STRUC) MUMPS_PAR
#else
      ! Fake declaration as DMUMPS_STRUC is shipped with MUMPS
      INTEGER MUMPS_PAR 
#endif
C
      character*1 anew_stif
#ifdef MUMPS5
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
C----------resol [K]{X}={F}------
      INTEGER  I  ,J ,ND,IP,ITMP,ICOV,ICALU,IER,ISETK0,IRIKS,
     .         F_DDL ,L_DDL
C     REAL
      my_real
     .  r2,u2,du2,tol,rr,ru,temp,fr,re,de,r1,nl_tol,r01,rr1,
     .  gapn,lamda,fa,uc2,um1,um2,pmax
C-----------------------------------------------
C------------------------------------------
c     insolv=1 =>  (elast) Newton modified  
c     insolv=2,3 =>  BFGS
c     insolv=4,5 =>  elastoplastic
C------------------------------------------
      fr = 0 ! was not always initialized before use
C------------------------------------------
C------------------------------------------
C     [K]:stiffness matrix [M]:preconditioning matrix
c        F_DDL=1+ITASK*NDDL/NTHREAD
c        L_ddl = (itask+1)*nddl/nthread
C
      tol =l_tol
      nl_tol=n_tol
C   
      lstol=ls_tol
      ip=iabs(inprint)
      de=0
      nd=3*numnod
      ru=ru0
      isetk0=isetk+icont
      IF (nitol==3.OR.nitol==13.OR.nitol==23
     .         .OR.nitol==123) THEN
       icalu=1
      ELSE
       icalu=0
      ENDIF
      iriks=0
      IF (ncycle>1.AND.idtc==3) THEN
       IF (ilast==0) iriks=1
      ENDIF
C----------------------
c      CALL MY_BARRIER
C---------------------
      isetk=0
      IF (icont>0) gapn=gap*zep9
      IF (imconv>=0) icont0=icont
C--------idtc=3 :modified Riks arc_length-----
      IF (it==0) THEN
       r01=sqrt(r02)
C--------ROLD=R2(IT=0)---------------------------------------
Ctmp       IF (ISPRB==1) R01=MIN(R01,SQRT(ROLD)/R01)
       re=one
C--------cas particulier---------------------------------------
       IF (r01*nl_tol<=em12.AND.nl_tol>em10.AND.idiv/=10) THEN
       IF (irefi/=5.OR.icont==0) THEN
        IF(inprint/=0)THEN
         WRITE(iout,*)
c         WRITE(IOUT,1008)IT,RE,R01,RE
         WRITE(iout,1108)it,anew_stif,re,r01,re
         WRITE(iout,*)
         IF(inprint<0)THEN
          WRITE(istdo,*)
c          WRITE(ISTDO,1008)IT,RE,R01,RE
          WRITE(istdo,1108)it,anew_stif,re,r01,re
          WRITE(istdo,*)
         ENDIF
        ENDIF
        relres(1) = re
        relres(2) = r01
        relres(3) = re
        idiv = -2       
        isetk=1
        RETURN
       ELSE
        f(1)=max(em20,f(1))
              r01=max(r01,em20)
       ENDIF
       ENDIF
       IF (insolv==4) icov = imconv
       imconv=0
       rr=rold
       ndiv=0
C---------
       IF (insolv==2.OR.insolv==3) CALL bfgs_0
C---------
       IF (ncycle==1) isign = 1
       IF (iriks>0) THEN
          CALL lin_solv(nddl  ,iddl  ,ndof  ,ikc   ,dg    ,
     1               dgr   ,tol   ,nnz   ,iadk  ,jdik  ,
     2               diag_k,lt_k  ,nddli ,iadi  ,jdii  ,   
     3               diag_i,lt_i  ,itok  ,iadm  ,jdim  ,   
     4               diag_m,lt_m  ,fext  ,de    ,inloc ,   
     5               fr_elem,iad_elem,w_ddl,itask0,icprec,
     6               istop ,a     ,ar     ,v    ,
     7               ms    ,x     ,ipari ,intbuf_tab ,
     8               num_imp,ns_imp,ne_imp,nsrem ,nsl  ,
     9               itc    ,graphe, itab, fac_k, ipiv_k,
     a               nk     ,nmonv,imonv  ,monvol,igrsurf,
     b               fr_mv ,volmon,ibfv  ,skew  ,
     c               xframe ,mumps_par,cddlp,ind_imp,rbid,
     e               irbe3 ,lrbe3 ,irbe2  ,lrbe2 )
 
C---------
         CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dg    ,dgr   ,u2    ,w_ddl  )
C-------distinquer converge or diverge---
citask0         IF (ITASK == 0) THEN
                IF (idiv/=-2) THEN
           CALL cp_real(nd,dd,dg0)
           IF (iroddl/=0) CALL cp_real(nd,ddr,dgr0) 
                ELSE
           dla_riks = isign*dt2
                ENDIF 
citask0         END IF !(ITASK == 0) THEN
          CALL produt_uhp2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                    dg    ,dgr   ,dg0   ,dgr0  ,uc2   ,
     .                    w_ddl )
C----------------------
c      CALL MY_BARRIER
C---------------------
citask0         IF (ITASK == 0) THEN
          CALL get_max(numnod,dg0,temp,i)
                um1=temp
          CALL get_max(numnod,dgr0,temp,i)
                um1=max(um1,temp,em20)
          CALL get_max(numnod,dg,temp,i)
                um2=temp
          CALL get_max(numnod,dgr,temp,i)
                um2=max(um2,temp,em20)
                uc2=uc2/um1/um2
C          LAMDA = ALEN/SQRT(U2)
               IF ((uc2+dla_riks/rfext)<zero) THEN
                isign = -1
               ELSE
                isign = 1
               ENDIF
c           print *,'UC2,DLA_RIKS,ISIGN=',UC2,DLA_RIKS,ISIGN
          IF (impdeb>0) THEN
           IF (ncycle>=ndeb0.AND.ncycle<=ndeb1) 
     .      WRITE(iout,1025)uc2,dla_riks,isign
                ENDIF
               IF (isign<0) THEN
                dla_riks = 2*isign*dt2
               ELSE
                dla_riks = zero
               ENDIF
C---------
citask0         END IF !(ITASK == 0) THEN
 
C----------------------
c      CALL MY_BARRIER
C---------------------
         lamda = isign*dt2
         IF (idiv/=-2.AND.icont==0)THEN
          CALL produt_vmhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                     dd    ,ddr   ,f     ,de    ,w_ddl )
          CALL integrator2_hp(iddl,ndof,ikc,d,dr,dd,ddr)
         ELSE 
          CALL frac_dd_hp(iddl,ndof,ikc,d,dr,dg,dgr,lamda)
          CALL produt_vmhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                     d     ,dr    ,f     ,de    ,w_ddl )
         ENDIF
       ELSE
       
        IF ((n_lim /= 1.OR.isprb /= 0).AND.
     .      (ncycle>1.AND.idiv/=-2.AND.icont==0.AND.ikt==0
     .     .OR.(ncycle==1.AND.idiv==10)))THEN
               IF (idtc==3.AND.ilast==1) THEN
          lamda = dla_riks
                dla_riks = zero
          CALL frac_d_hp(iddl,ndof,ikc,dd,ddr,lamda)
               ENDIF
         CALL produt_vmhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                   dd    ,ddr   ,f     ,de    ,w_ddl )
        ELSE
         IF (insolv==4.AND.icov>=0) imconv=-1
 
         CALL lin_solv(nddl  ,iddl  ,ndof  ,ikc   ,dd    ,
     1               ddr   ,tol   ,nnz   ,iadk  ,jdik  ,
     2               diag_k,lt_k  ,nddli ,iadi  ,jdii  ,   
     3               diag_i,lt_i  ,itok  ,iadm  ,jdim  ,   
     4               diag_m,lt_m  ,f     ,de    ,inloc ,   
     5               fr_elem,iad_elem,w_ddl,itask0,icprec,
     6               istop ,a     ,ar     ,v    ,
     7               ms    ,x     ,ipari ,intbuf_tab ,
     8               num_imp,ns_imp,ne_imp,nsrem ,nsl  ,
     9               itc    ,graphe, itab, fac_k, ipiv_k,
     a               nk     ,nmonv,imonv  ,monvol,igrsurf,
     b               fr_mv ,volmon,ibfv  ,skew  ,
     c               xframe ,mumps_par,cddlp,ind_imp,rbid,
     e               irbe3 ,lrbe3 ,irbe2  ,lrbe2 )
         IF (de<zero.AND.irefi>1.AND.irefi<5) imconv=-2
         IF (icont>0.AND.idsgap>0) THEN
         CALL get_max(numnod,dd,pmax,i)
C----------------------
c      CALL MY_BARRIER
C---------------------
          IF (pmax>gapn) THEN
           temp=gapn/pmax
           CALL frac_d_hp(iddl,ndof,ikc,dd,ddr,temp)
                 de = de * temp
          ENDIF
         ENDIF
        ENDIF
C
        CALL integrator2_hp(iddl,ndof,ikc,d,dr,dd,ddr)
       ENDIF !IF (IRIKS>0)
C----------------------
c      CALL MY_BARRIER
C---------------------
citask0       IF (ITASK == 0) THEN
       IF (imconv>=0) THEN
        CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  d    ,dr    ,u2    ,w_ddl )
        eimp=eimp+de
        IF (de>ep10.AND.iline==0.AND.idyna==0
     .     .AND.ncycle==1) THEN
         IF(ispmd==0)THEN
          WRITE(iout,1030)eimp  
          WRITE(istdo,1030)eimp  
         ENDIF
               CALL imp_stop(-1)
        ENDIF
C--------for idtc>=2-----
        u02=u2
        IF(inprint/=0)THEN
         ru=sqrt(u2)
c         WRITE(IOUT,*)  
c         WRITE(IOUT,1002)IT,RU,R01,EIMP
         WRITE(iout,1102)
         WRITE(iout,1003)it,anew_stif,ru,r01,eimp
          IF(inprint<0)THEN
c          WRITE(ISTDO,*)  
c          WRITE(ISTDO,1002)IT,RU,R01,EIMP
          WRITE(istdo,1102)
          WRITE(istdo,1003)it,anew_stif,ru,r01,eimp
         ENDIF
        ENDIF
 
        relres(1) = ru
        relres(2) = r01
        relres(3) = eimp
       ELSE
          WRITE(iout,*)  
          WRITE(iout,1013) de
          IF(inprint<0)THEN
           WRITE(istdo,*)  
           WRITE(istdo,1013)de
          ENDIF
C
       ENDIF !IF (IMCONV>=0)
citask0       END IF !(ITASK == 0) THEN
C---------IT>0---------
      ELSE
C--------R=Fext-Fint correction due to Fext---------
       IF (iriks>0.AND.dla_riks/=zero) THEN
        CALL vaxpy_hp(nddl, f   ,fext,dla_riks)
       ENDIF
       CALL produt_hp(nddl,f,f,w_ddl,r2)                      
C----------------------
c      CALL MY_BARRIER
C---------------------
      IF (r2>=zero.AND.r2<ep30) THEN
       rr=sqrt(r2/r02)
       IF (it==1.AND.isprb==1.AND.idiv/=-2) THEN
        rr1=sqrt(r2/rold)
C
        IF (icont>0.AND.rr<one) rr1=rr
       ELSE
        rr1=rr
       ENDIF
       re=de0/eimp
       IF (it==1) THEN
        ru=one
       ELSEIF (icalu>0) THEN
         CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                   d     ,dr    ,u2    ,w_ddl )
         CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                   dd    ,ddr   ,du2   ,w_ddl ) 
C----------------------
c      CALL MY_BARRIER
C---------------------
         ru=ls(3)*sqrt(du2/u2)
         IF (insolv==5) ru=sqrt(du2/u2)
       ENDIF 
       CALL crit_ite(it,ru,rr1,re,ndiv,nl_tol)
      ELSE
       imconv =-2  
       rr1 =ep30
      ENDIF
C----------------------
c      CALL MY_BARRIER
C---------------------
       IF (imconv==-3) THEN
        IF (ncycle==1.OR.idiv==-2) imconv=-2
        r02=r2
       ENDIF
C----------------------
c      CALL MY_BARRIER
C---------------------
       IF (imconv==0) THEN
         idiv=0
C-------start line-search------------
        IF (iline_s/=1) THEN
         ls(1)=one
         ls(2)=zero
        ENDIF
         ls(3)=one
         IF (ndiver>0) r02=r2
       ENDIF
       IF(imconv==1) THEN
        istop=0
        idiv=0
C------For time correction-with Riks----   
         IF (it>1.AND.icalu==0) THEN
          CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                    d     ,dr    ,u2    ,w_ddl )
          CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                    dd    ,ddr   ,du2   ,w_ddl )
C----------------------
c      CALL MY_BARRIER
C---------------------
          ru=ls(3)*sqrt(du2/u2)
          u02=u2
         ELSEIF (idtc>=2) THEN
          CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                    d     ,dr    ,u02   ,w_ddl )
C----------------------
c      CALL MY_BARRIER
C---------------------
         ENDIF
C
        IF(inprint/=0)THEN
c         WRITE(IOUT,*)
c         WRITE(IOUT,1008)IT,RU,RR,RE
         WRITE(iout,1108)it,anew_stif,ru,rr,re
         WRITE(iout,*)
         IF(inprint<0)THEN
c          WRITE(ISTDO,*)
c          WRITE(ISTDO,1008)IT,RU,RR,RE
          WRITE(istdo,1108)it,anew_stif,ru,rr,re
          WRITE(istdo,*)
         ENDIF
        ENDIF
        relres(1) = ru
        relres(2) = rr
        relres(3) = re
       ELSEIF(imconv<=-2) THEN
C--------diverge----
C        write(iout,*)'diverge',IMCONV,IT ,IDIV  
        IF(inprint/=0)THEN
         WRITE(iout,*)
         IF (rr1>one) THEN
          WRITE(iout,1011)rr1
         ELSEIF(it>nl_dtn) THEN
          WRITE(iout,1010)
         ELSEIF(idtc==3) THEN
          WRITE(iout,1020)
         ENDIF
         WRITE(iout,*)
         IF(inprint<0)THEN
          WRITE(istdo,*)
          IF (rr1>one) THEN
           WRITE(istdo,1011)rr1
          ELSEIF(it>nl_dtn) THEN
           WRITE(istdo,1010)
          ELSEIF(idtc==3) THEN
           WRITE(istdo,1020)
          ENDIF
          WRITE(istdo,*)
         ENDIF
         IF(imconv==-2) THEN
          WRITE(iout,1012)
          IF(inprint<0)WRITE(istdo,1012)
         ENDIF
        ENDIF
        relres(1) = ru
        relres(2) = rr1
        relres(3) = re
       ELSE
C--------CAS IMCONV=0,-1---add line-search for constant dt-
c        IF (IDTC>0) THEN
        IF (iline_s>0.AND.isign>=0.AND.irwall==0) THEN
         fr=ls(3)
         IF (nitol/=2.AND.nitol/=4.OR.iline_s==1) THEN
C----------DD,DDR are modified in some cases         
          IF (ls1(3)/=zero ) THEN
           de = ls1(3)
          ELSE
          CALL produt_vmhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                     dd    ,ddr   ,f     ,de    ,w_ddl  )
          END IF
C----------------------
c      CALL MY_BARRIER
C---------------------
          r1=de/de0
          IF (iline_s==3) THEN
           r1=abs(de/de0)
           IF (rr>one) r1=one
           IF (rr>one.AND.irefi>=2) r1=ls(1)
           temp=ep02
          ENDIF
         ELSE
          r1=rr/rold
          temp=ep03
         ENDIF
         icov = imconv
C------------
citask0        IF (ITASK == 0) THEN
C------------
         IF (iline_s==3) THEN
          CALL line_s(ls(1),ls(2),r1,fr,lstol,idiv,icont,temp)
         ELSEIF (iline_s==1) THEN
          CALL line_s1(r1,fr,ls(1),ls(2),ls1(1),ls1(2),idiv,lstol,
     .                 icont,icont0,iriks)
         ELSEIF (iline_s==2) THEN
          CALL line_s(ls(1),ls(2),r1,fr,lstol,idiv,icont,temp)
         ENDIF
         IF (impdeb>0.AND.ispmd==0) THEN
          IF (ncycle>=ndeb0.AND.ncycle<=ndeb1) then
           WRITE(iout,*)'R1,FR,IDIV=',r1,fr,idiv
C           WRITE(IOUT,*)'DE,DE0=',DE,DE0
          end if
         ENDIF
C------------
citask0        END IF !(ITASK == 0) THEN
C----------------------
c       CALL MY_BARRIER
C---------------------
C
         IF ((insolv==2.OR.insolv==3).AND.fr>one.AND.ikt==0) 
     .       fr = min(fr,two)
        ENDIF
C        IF (IABS(INSOLV)==5) THEN
C----------suprimed----------
C        ELSE  
        IF (imconv==0) THEN
C----------schema normale---------
         IF (it>1) THEN
C------------
          CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                    d    ,dr    ,u2    ,w_ddl  )
          CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                    dd    ,ddr   ,du2   ,w_ddl  )
C----------------------
c      CALL MY_BARRIER
C---------------------
          ru=ls(3)*sqrt(du2/u2)
         ENDIF
C------------
citask0        IF (ITASK == 0) THEN
C------------
         IF (insolv==2.OR.insolv==3) CALL bfgs_ls(ls(3))
         IF (rr>=rold) THEN
          ndiv = ndiv +1
         ELSE
          ndiv = 0
         ENDIF
         IF (inprint/=0)THEN
          IF(mod(it,ip)==0)THEN
c           WRITE(IOUT,*)  
           WRITE(iout,1003)it,anew_stif,ru,rr,re
           IF(inprint<0)THEN
c            WRITE(ISTDO,*)  
            WRITE(istdo,1003)it,anew_stif,ru,rr,re
           ENDIF
          ENDIF
         ENDIF
         relres(1) = ru
         relres(2) = rr
         relres(3) = re
         IF(itc>=n_lim.AND.ismdisp==0) THEN
c          IF(INPRINT/=0)THEN
c            WRITE(IOUT,*)  
c            WRITE(IOUT,1006)ITC
c            IF(INPRINT<0)THEN
c             WRITE(ISTDO,*)  
c             WRITE(ISTDO,1006)ITC
c            ENDIF
c          ENDIF
          itc=0
         ENDIF
C------------
citask0        END IF !(ITASK == 0) THEN
C----------------------
c      CALL MY_BARRIER
C---------------------
         IF (iriks>0) THEN
          IF (isetk0==1) THEN
           CALL lin_solv(nddl  ,iddl  ,ndof  ,ikc   ,dg    ,
     1               dgr   ,tol   ,nnz   ,iadk  ,jdik  ,
     2               diag_k,lt_k  ,nddli ,iadi  ,jdii  ,   
     3               diag_i,lt_i  ,itok  ,iadm  ,jdim  ,   
     4               diag_m,lt_m  ,fext  ,de    ,inloc ,   
     5               fr_elem,iad_elem,w_ddl,itask0,icprec,
     6               istop ,a     ,ar     ,v    ,
     7               ms    ,x     ,ipari ,intbuf_tab ,
     8               num_imp,ns_imp,ne_imp,nsrem ,nsl  ,
     9               itc    ,graphe, itab, fac_k, ipiv_k,
     a               nk     ,nmonv,imonv  ,monvol,igrsurf,
     b               fr_mv ,volmon,ibfv  ,skew  ,
     c               xframe ,mumps_par,cddlp,ind_imp,rbid,
     e               irbe3 ,lrbe3 ,irbe2  ,lrbe2 )
                icprec = 0
                idsc = 0
          END IF 
          CALL produt_vmhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                     dg    ,dgr   ,f     ,um1   ,w_ddl )
         END IF 
C
        CALL lin_solv(nddl  ,iddl  ,ndof  ,ikc   ,dd    ,
     1                 ddr   ,tol   ,nnz   ,iadk  ,jdik  ,
     2                 diag_k,lt_k  ,nddli ,iadi  ,jdii  ,   
     3                 diag_i,lt_i  ,itok  ,iadm  ,jdim  ,   
     4                 diag_m,lt_m  ,f     ,de    ,inloc ,   
     5                 fr_elem,iad_elem,w_ddl,itask0,icprec,
     6                 istop ,a     ,ar     ,v    ,
     7                 ms    ,x     ,ipari ,intbuf_tab ,
     8                 num_imp,ns_imp,ne_imp,nsrem ,nsl  ,
     9                 itc    ,graphe, itab, fac_k, ipiv_k,
     a                 nk     ,nmonv,imonv  ,monvol,igrsurf,
     b                 fr_mv ,volmon,ibfv  ,skew  ,
     c                 xframe ,mumps_par,cddlp,ind_imp,rbid,
     d                 irbe3 ,lrbe3 ,irbe2  ,lrbe2 )
         IF (icont>0.AND.idsgap>0) THEN
          CALL get_max(numnod,dd,pmax,i)
C----------------------
c      CALL MY_BARRIER
C---------------------
          IF (pmax>gapn) THEN
           temp=gapn/pmax
           CALL frac_d_hp(iddl,ndof,ikc,dd,ddr,temp)
                 de = de * temp
          ENDIF
         ENDIF
C----------------------
c      CALL MY_BARRIER
C---------------------
C---------Riks type arc length------   
         IF (iriks>0) THEN
                lamda = dla_riks
                IF (ial_m==1) THEN
            CALL al_constraint1_hp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                          dd    ,ddr   ,dg    ,dgr   ,d     ,
     .                          dr    ,w_ddl ,alen  ,lamda ,scal_riks,
     .                          ier   )
C----------------------
c      CALL MY_BARRIER
C---------------------
          IF (ier==1) THEN
                 imconv=-2
            IF(inprint/=0)THEN
             WRITE(iout,*)
             WRITE(iout,1020)
             WRITE(iout,*)
             IF(inprint<0)THEN
              WRITE(istdo,*)
              WRITE(istdo,1020)
              WRITE(istdo,*)
             ENDIF
            ENDIF
           END IF !(IER==1) THEN
          ELSEIF (ial_m==2) THEN
            CALL al_constraint2_hp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                          dd    ,ddr   ,dg    ,dgr   ,d     ,
     .                          dr    ,w_ddl ,alen  ,lamda ,scal_riks)
C----------------------
c      CALL MY_BARRIER
C---------------------
            END IF
            IF (imconv>=0) THEN
           CALL frac_dd_hp(iddl,ndof,ikc,d ,dr ,dg,dgr,lamda)
citask0           IF (ITASK == 0) THEN
                  dla_riks = dla_riks + lamda
            eimp=eimp+dla_riks*um1
citask0           END IF !(ITASK == 0) THEN
C----------------------
c      CALL MY_BARRIER
C---------------------
          ENDIF
         END IF !IF (IRIKS>0) THEN
C   
         CALL integrator2_hp(iddl,ndof,ikc,d,dr,dd,ddr)
         eimp=eimp+de
c   if (de<0) print *,'****negative de=',de
C
        ELSEIF(imconv==-1) THEN
C----------line-search---------
         temp=-ls(3)+fr
        CALL frac_dd_hp(iddl,ndof,ikc,d,dr,dd,ddr,temp)
C----------------------
c      CALL MY_BARRIER
C---------------------
        ls(3) = fr
        END IF !IF (IMCONV==0) THEN
C      
       END IF !IF(IMCONV==1) THEN
C-----fin if IT=0
      ENDIF
      IF (imconv>=0) THEN
       ru0=ru
       de0=de
       IF (abs(de0)<em20) de0=em20
       rold=max(em20,rr)
      ENDIF
C-----REFORMING [K]----      
      IF (ismdisp>0) THEN
       IF (itc>=n_lim.AND.imconv>=0) THEN
         itc = 0
         isetk=1
       END IF
       IF (itc==0.AND.imconv>=0) THEN
        IF (ncycle==1.OR.n_lim==1) isetk=1
        IF (isolv==5.OR.isolv==6) isetk=1
       END IF
       IF (irwall >0 .AND.imconv == 1) isetk=1
       IF ((itc==0.OR.imconv==1).AND.ikt>0) isetk=1
      ELSE
C----condition IMCONV==1 reforming [K] for nothing (no resolution)-       
       IF (itc==0.OR.imconv==1) isetk=1
      END IF
      IF ((idtc==3.OR.ikt>0).AND.imconv<=-2) isetk=1
C-----if still diverge, try w/o [K] reforming      
      IF (imconv<=-2.AND.rr1>ep20.AND.idiv0>=0) isetk=1
      IF ((istop==1.OR.istop==2).AND.(isolv==5.OR.isolv==6))THEN 
        istop = 0
        isetk = 1
        imconv=-2
      ENDIF
      IF (imconv<=-2.AND.isprb>0.AND.rr1>ep04) isetk=1
! sb 
      IF(imconv<= -2.AND.n_lim == 1 .AND. isprb == 0)  isetk = 1
      
C----------------------
c      CALL MY_BARRIER
C---------------------
Ctmp      IF (ITC==0.OR.IMCONV==1.OR.IMCONV<=-2) ISETK=1
c 1002 FORMAT(5X,'NL_ITERATION=',I5,1X,'  INITIAL RESIDUAL NORM=',3E11.4)
 1102 FORMAT(3x,78('-')/
     .       12x,'Stif. Mat.'/
     .       6x,'Iter',2x, ' reformed ',5x,'|du|/|u|',3x,'|r|/|r0|',3x,'|dE|/|E|',1x,'Conv.stat.'/
     .       3x,78('-'))
c 1003 FORMAT(5X,'NL_ITERATION=',I5,1X,' RELATIVE RESIDUAL NORM=',3E11.4)
 1003 FORMAT(5x,i5,6x,a1,5x,3(1x,1pe10.3))
 1005 FORMAT(3x,'TOLERANCE FOR LINEAR ITERATIVE SOLVER :',2x,e11.4)
 1006 FORMAT(3x,'--STIFFNESS MATRIX IS RESET AFTER ',i4,
     .          ' ITERATIONS--')
 1007 FORMAT(3x,'--ITERATION DIVERGE with R=',e11.4,2x,
     .       'STIFFNESS MATRIX WILL BE RESET')
c 1008 FORMAT(3X,'* CONVERGED WITH ',I4,1X,'ITERATIONS,',
c     .          ' |du|/|u|,|r|/|r0|,|dE|/|E|=',3E11.4)
 1108 FORMAT(5x,i5,6x,a1,5x,3(1x,1pe10.3),5x, 'C')
 1009 FORMAT(3x,'--ITERATION DIVERGE with RELATIVE R=',e11.4/,
     .       3x,'RESET ITERATION WITH STABILIZATION BY SCALING= ',e11.4)
 1010 FORMAT(3x,'--ITERATION DIVERGE with MAX_ITER REACHED--')
 1011 FORMAT(3x,'ITERATION DIVERGE with RELATIVE R=',e11.4)
 1012 FORMAT(3x,'--RESET ITERATION WITH NEW TIMESTEP--')
 1013 FORMAT(3x,'ITERATION DIVERGE with NEGATIVE ENERGY DE=',e11.4)
 1020 FORMAT(3x,'--ITERATION DIVERGE with RIKS Method--')
 1025 FORMAT(3x,'UC2,DLA_RIKS,ISIGN=',2e11.4,i4)
 1030 FORMAT(3x,'CHECK CONSTRAINT CONDITIONS, TOO LARGE ENERGY VALUE=',
     .           2x,e11.4)
      RETURN
#endif
      END
 
!||====================================================================
!||    crit_ite   ../engine/source/implicit/nl_solv.F
!||--- called by ------------------------------------------------------
!||    nl_solv    ../engine/source/implicit/nl_solv.F
!||====================================================================
      SUBROUTINE crit_ite(IT,UR,RR,ER,NDIV,TOL)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "impl1_c.inc" 
#include      "impl2_c.inc" 
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  IT,NDIV
C     REAL
      my_real
     .  UR,RR,ER,TOL
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      my_real
     .  err,tolr,told
C--------cas particulier---------------------------------------
C      TOLD=EP04
C      IF (ILINE_S==1) TOLD=EP03 
C      IF (ISMDISP==1) TOLD=EP10 
      told = tol_div
      tolr = one+tol
C--------1:E,2:r,3:E+U,4:r+U,5,r+E,6,r+E+U-----
C--new:  1:E,2:r,3:U,12:E+r,13,E+U,23,r+U,123,r+E+U
      IF (nitol==1) THEN
       err=abs(er)/tol
       IF (it==1.OR.rr>=one) err=err+tolr 
      ELSEIF (nitol==2) THEN
       err=rr/tol
      ELSEIF (nitol==3) THEN
       err=ur/tol
C       ERR=MAX(SQRT(ABS(ER)),UR)
       IF (rr>=one) err=err+tolr 
      ELSEIF (nitol==12) THEN
       err=max(rr/n_tolf,abs(er)/n_tole)
c       ERR=MAX(RR,UR)
      ELSEIF (nitol==13) THEN
       err=max(ur/n_tolu,abs(er)/n_tole)
c       ERR=MAX(SQRT(ABS(ER)),RR)
      ELSEIF (nitol==23) THEN
       err=max(ur/n_tolu,rr/n_tolf)
c       ERR=MAX(SQRT(ABS(ER)),RR,UR)
      ELSEIF (nitol==123) THEN
       err=max(ur/n_tolu,rr/n_tolf,abs(er)/n_tole)
      ENDIF
C
       IF (err<=one) THEN
C       IF (ERR<=TOL) THEN
        imconv=1
!sb       ELSEIF (UR<MIN(EM05,TOL).AND.RR<EM01) THEN
       ELSEIF ((n_lim /= 1.OR.isprb /= 0).AND.ur<min(em05,tol).AND.rr<em01) THEN
       imconv=1
       ELSEIF (imconv<0) THEN
C--------line-search-----
       ELSE
        imconv=0
C------divergence-----------
!sb
       IF(n_lim /= 1 .OR. isprb /= 0) THEN
        IF (rr>tolr) THEN
         IF (it==1) THEN
          imconv=-3
          IF (ndiv<ndiver) imconv=0
         ELSE
          imconv=-2
          IF ((ndiv<ndiver.OR.ilast==1).AND.rr<told) imconv=0
         ENDIF
        ENDIF
       ENDIF
 
       IF (idtc>0.AND.it>nl_dtn) imconv=-2
       ENDIF
C-----------------------------------------------
      RETURN
      END
!||====================================================================
!||    line_s    ../engine/source/implicit/nl_solv.F
!||--- called by ------------------------------------------------------
!||    nl_solv   ../engine/source/implicit/nl_solv.F
!||====================================================================
      SUBROUTINE line_s(R0,S0,R1,S,DTOL,IDIV,IINT,PREC)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "impl1_c.inc" 
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER IDIV,IINT
C     REAL
      my_real
     .  r0,r1,s0,s,dtol,prec
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .  s1,dr,r
C-----------------------------------------------
C      R1=ABS(E1/E0)
C      IF (R1<=ONE) THEN
      dr=abs(r0-r1)
      r=min(r1/prec,dr)
c      R=MIN(R1,DR*PREC)
      IF ((r<dtol.OR.idiv<-nls_lim) 
     .    .OR.(r1>one.AND.iint>0.AND.irefi<2)) THEN
       IF (imconv==-1)imconv=0 
c       IDIV=0
       idiv=idiv-1
       RETURN
      ELSEIF (r1>r0.AND.imconv==-1) THEN
       idiv=-nls_lim-2
       s=s0
       RETURN
      ENDIF
      imconv=-1     
      s1=s
      s=abs(s0*r1-s1*r0)/dr
      IF (r1>one) THEN
C----------diverge---------------------------
       IF (s>one) s=half*s1/r1
       s=max(em01,s)
C       IDIV=IDIV+1
      ELSE
C----- extrapolation --------
       IF (s>max(s1,s0).OR.s<min(s1,s0)) THEN
        s=min(three_half*s1,s)
        s=min(five,s)
       ENDIF
      ENDIF
      s0=s1
      r0=r1
      idiv=idiv-1
C
      RETURN
      END
!||====================================================================
!||    get_max      ../engine/source/implicit/nl_solv.F
!||--- called by ------------------------------------------------------
!||    nl_solv      ../engine/source/implicit/nl_solv.F
!||--- calls      -----------------------------------------------------
!||    spmd_max_s   ../engine/source/mpi/implicit/imp_spmd.F
!||====================================================================
      SUBROUTINE get_max(N,D,VMAX,I)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER N,I
      my_real
     .  D(3,*),VMAX
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER K
      my_real
     .  t
C-----------------------------------------------
       vmax=zero
       i=1
       DO k=1,n
        t=max(abs(d(1,k)),abs(d(2,k)),abs(d(3,k)))
        IF (t>vmax) THEN
         vmax = t
         i=k
        ENDIF
       ENDDO 
       IF (nspmd>1) CALL spmd_max_s(vmax)
C
      RETURN
      END
!||====================================================================
!||    line_s1   ../engine/source/implicit/nl_solv.F
!||--- called by ------------------------------------------------------
!||    nl_solv   ../engine/source/implicit/nl_solv.F
!||====================================================================
      SUBROUTINE line_s1(E1,S,EP,SP,EN,SN,IDIV,DTOL,ICONT,ICONT0,IRIKS)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "impl1_c.inc" 
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER IDIV,ICONT,ICONT0,IRIKS
C     REAL
      my_real
     .  e1,s,ep,sp,en,sn,dtol
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,IM
      my_real
     .  sold,stmp,s2,tol,eold
C-------min e=r.du-----------------E1-> already relative /E0-----------------------
C------find sl so that el<0 
       imconv =-1
       eold=zero
       im = max(icont,icont0)
       IF (e1>zero) THEN
        IF (idiv==0) THEN
         idiv = idiv + 1
         sp = one
         ep = e1
         s = s + one
         IF (iriks>0.OR.im > 0) THEN
          s = one
          imconv =0
         END IF
         RETURN
        ELSEIF (idiv>0) THEN
         sold = s
         IF (abs(e1-ep)<dtol.OR.e1>ep.OR.idiv>nls_lim) THEN   
          IF (e1>ep) s = s - one
          imconv =0
         ELSE
          s = s + one
         ENDIF
         idiv = idiv + 1
         sp = sold
         ep = e1
         RETURN
        ENDIF
        sp = s
        ep = e1
       ELSE
        IF (idiv==0) THEN   
         sp = zero
         ep = one
         s = half
         idiv = idiv - 1
         sn = one
         en = e1
         RETURN
        ELSEIF (sp==zero) THEN
          IF (abs(e1-en)<dtol.OR.idiv<-nls_lim) THEN   
           imconv =0
          ELSEIF (iriks>0.OR.ismdisp>0) THEN
           imconv =0
          ELSE
c          SN = S
c    EN = E1
           s = s*half
          ENDIF
          idiv = idiv - 1
          RETURN
        ELSEIF (idiv>0) THEN
         idiv = -idiv
        ENDIF
        eold=en
        sn = s
        en = e1
       ENDIF
       idiv = idiv - 1
c----------return to       
       sold = s
       im = max(icont,icont0)
       IF (im>0) THEN
c        S2 = (EP-EN)/SN+SN
c        STMP = HALF*(S2-SQRT(S2*S2-FOUR*EP))
        stmp = half*(-sp+sn)
       ELSE
        stmp = ep*(sn-sp)/(ep-en)
       ENDIF
       s = sp + stmp
       stmp =(sold-s)/sold
       stmp = max(abs(stmp),abs(e1-eold))
       IF (stmp<dtol.OR.idiv<-nls_lim) THEN
        imconv =0
       ELSEIF (iriks>0.AND.
     .         (s>two.OR.s<em01)) THEN
        imconv =0
       ELSE
        imconv =-1
       ENDIF
C
      RETURN
      END
!||====================================================================
!||    al_constraint1   ../engine/source/implicit/nl_solv.F
!||--- calls      -----------------------------------------------------
!||    frac_dd          ../engine/source/implicit/integrator.F
!||    integrator2      ../engine/source/implicit/integrator.F
!||    produt_u         ../engine/source/implicit/produt_v.F
!||    produt_u2        ../engine/source/implicit/produt_v.F
!||====================================================================
      SUBROUTINE al_constraint1(NDDL0 ,NDDL  ,IDDL  ,NDOF  ,IKC   ,
     .                          DD    ,DDR   ,DG    ,DGR   ,DI    ,
     .                          DIR   ,W_DDL ,L_A   ,LAMDA ,SW2   ,
     .                          IER   )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com04_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  NDDL,NDDL0,IDDL(*)  ,NDOF(*)  ,IKC(*),W_DDL(*) ,IER
C     REAL
      my_real
     .  LAMDA,DG(3,*),DGR(3,*),DI(3,*),DIR(3,*),DD(3,*),DDR(3,*),
     .  L_A,SW2
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .  A,B,C,UI2,UIUD,UIUT,UT2,LAM1,LAM2,S,TH1,TH2,FAC
C-----------------------------------------------
         IER=0
         call produt_u(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                 dg    ,dgr   ,ut2   ,w_ddl )
         a = ut2+sw2
C-------DD< DI+DD   
         CALL integrator2(1,numnod,iddl,ndof,ikc,dd,ddr,di,dir)
         CALL produt_u2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dg    ,dgr   ,dd    ,ddr   ,uiut  ,
     .                  w_ddl )
         b= two*(uiut+sw2*lamda)
         CALL produt_u(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                 dd    ,ddr   ,ui2   ,w_ddl )
         c = ui2-l_a*l_a+sw2*lamda
           s= b*b-four*a*c
          IF (s>=0) THEN
           s=sqrt(s)
           lam1=half*(-b+s)/a
           lam2=half*(-b-s)/a
         CALL produt_u2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dd    ,ddr   ,di    ,dir   ,uiud  ,
     .                  w_ddl )
         CALL produt_u2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dg    ,dgr   ,di    ,dir   ,uiut  ,
     .                  w_ddl )
         th1= uiud+lam1*uiut
         th2= uiud+lam2*uiut
           IF (th1>zero.AND.th2<zero) THEN
            lamda =lam1
           ELSEIF (th1<zero.AND.th2>zero) THEN
            lamda =lam2
           ELSEIF (th1>zero.AND.th2>zero) THEN
            lamda =min(lam1,lam2)
           ELSE
            ier=1
           ENDIF
        ELSE
           ier=1
          ENDIF
           fac =-one
        CALL frac_dd(1,numnod,iddl,ndof,ikc,dd,ddr,di,dir,fac)
C
      RETURN
      END
!||====================================================================
!||    al_constraint2   ../engine/source/implicit/nl_solv.F
!||--- calls      -----------------------------------------------------
!||    produt_u         ../engine/source/implicit/produt_v.F
!||    produt_u2        ../engine/source/implicit/produt_v.F
!||====================================================================
      SUBROUTINE al_constraint2(NDDL0 ,NDDL  ,IDDL  ,NDOF  ,IKC   ,
     .                          DD    ,DDR   ,DG    ,DGR   ,DI    ,
     .                          DIR   ,W_DDL ,L_A   ,LAMDA ,SW2   )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  NDDL,NDDL0,IDDL(*)  ,NDOF(*)  ,IKC(*),W_DDL(*)  
      my_real
     .  lamda,dg(3,*),dgr(3,*),di(3,*),dir(3,*),dd(3,*),ddr(3,*),
     .  l_a,sw2
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .  gi,ui2,uiud,uiut,s
C-----------------------------------------------
         CALL produt_u(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                 di    ,dir   ,ui2   ,w_ddl )
         gi = sqrt(ui2)
         CALL produt_u2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dg    ,dgr   ,di    ,dir   ,uiut  ,
     .                  w_ddl )
         CALL produt_u2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dd    ,ddr   ,di    ,dir   ,uiud  ,
     .                  w_ddl )
         s = sw2*lamda   
         lamda = ((l_a-gi)*gi-uiud)/(uiut+s)
C
      RETURN
      END
!||====================================================================
!||    b_correct   ../engine/source/implicit/nl_solv.F
!||====================================================================
      SUBROUTINE b_correct(NDDL  ,F   ,FEXT  ,DLAMDA )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  NDDL 
      my_real
     .  dlamda,f(*),fext(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
C-----------------------------------------------
         DO I =1,nddl
            f(i)=f(i)+dlamda*fext(i)
           ENDDO
C
      RETURN
      END
!||====================================================================
!||    b_correct_h   ../engine/source/implicit/nl_solv.F
!||====================================================================
      SUBROUTINE b_correct_h(F_DDL  ,L_DDL  ,F   ,FEXT  ,DLAMDA )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  F_DDL  ,L_DDL 
C     REAL
      my_real
     .  dlamda,f(*),fext(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
C-----------------------------------------------
         DO i =f_ddl  ,l_ddl
                f(i)=f(i)+dlamda*fext(i)
               ENDDO
C
      RETURN
      END
!||====================================================================
!||    al_constrainth1   ../engine/source/implicit/nl_solv.F
!||--- calls      -----------------------------------------------------
!||    frac_dd           ../engine/source/implicit/integrator.F
!||    integrator2       ../engine/source/implicit/integrator.F
!||    my_barrier        ../engine/source/system/machine.F
!||    produt_uh         ../engine/source/implicit/produt_v.F
!||    produt_uh2        ../engine/source/implicit/produt_v.F
!||====================================================================
      SUBROUTINE al_constrainth1(NDDL0 ,NDDL  ,IDDL  ,NDOF  ,IKC   ,
     .                          DD    ,DDR   ,DG    ,DGR   ,DI    ,
     .                          DIR   ,W_DDL ,L_A   ,LAMDA ,SW2   ,
     .                          FT_N  ,LT_N  ,F_DDL ,L_DDL ,ITASK ,
     .                          IER   )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  NDDL,NDDL0,IDDL(*)  ,NDOF(*)  ,IKC(*),W_DDL(*) ,IER,
     .         FT_N  ,LT_N  ,F_DDL ,L_DDL ,ITASK
C     REAL
      my_real
     .  lamda,dg(3,*),dgr(3,*),di(3,*),dir(3,*),dd(3,*),ddr(3,*),
     .  l_a,sw2
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .  A,B,C,UI2,UIUD,UIUT,UT2,LAM1,LAM2,S,TH1,TH2,FAC
C-----------------------------------------------
         IER=0
         call produt_uh(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                 dg    ,dgr   ,ut2   ,w_ddl  ,f_ddl ,
     .                 l_ddl ,itask )
C----------------------
      CALL my_barrier
C---------------------
         a = ut2+sw2
C-------DD< DI+DD   
         CALL integrator2(ft_n,lt_n,iddl,ndof,ikc,dd,ddr,di,dir)
C----------------------
      CALL my_barrier
C---------------------
         CALL produt_uh2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dg    ,dgr   ,dd    ,ddr   ,uiut  ,
     .                  w_ddl ,f_ddl ,l_ddl ,itask )
         CALL produt_uh(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                 dd    ,ddr   ,ui2   ,w_ddl  ,f_ddl ,
     .                 l_ddl ,itask )
C----------------------
      CALL my_barrier
C---------------------
         b= two*(uiut+sw2*lamda)
         c = ui2-l_a*l_a+sw2*lamda
           s= b*b-four*a*c
          IF (s>=0) THEN
           s=sqrt(s)
           lam1=half*(-b+s)/a
           lam2=half*(-b-s)/a
         CALL produt_uh2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dd    ,ddr   ,di    ,dir   ,uiud  ,
     .                  w_ddl ,f_ddl ,l_ddl ,itask )
         CALL produt_uh2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dg    ,dgr   ,di    ,dir   ,uiut  ,
     .                  w_ddl ,f_ddl ,l_ddl ,itask )
C----------------------
      CALL my_barrier
C---------------------
         th1= uiud+lam1*uiut
         th2= uiud+lam2*uiut
           IF (th1>zero.AND.th2<zero) THEN
            lamda =lam1
           ELSEIF (th1<zero.AND.th2>zero) THEN
            lamda =lam2
           ELSEIF (th1>zero.AND.th2>zero) THEN
            lamda =min(lam1,lam2)
           ELSE
            ier=1
           ENDIF
        ELSE
           ier=1
          ENDIF
           fac =-one
        CALL frac_dd(ft_n,lt_n,iddl,ndof,ikc,dd,ddr,di,dir,fac)
C
      RETURN
      END
!||====================================================================
!||    al_constrainth2   ../engine/source/implicit/nl_solv.F
!||--- calls      -----------------------------------------------------
!||    my_barrier        ../engine/source/system/machine.F
!||    produt_uh         ../engine/source/implicit/produt_v.F
!||    produt_uh2        ../engine/source/implicit/produt_v.F
!||====================================================================
      SUBROUTINE al_constrainth2(NDDL0 ,NDDL  ,IDDL  ,NDOF  ,IKC   ,
     .                          DD    ,DDR   ,DG    ,DGR   ,DI    ,
     .                          DIR   ,W_DDL ,L_A   ,LAMDA ,SW2   ,
     .                          FT_N  ,LT_N  ,F_DDL ,L_DDL ,ITASK )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  NDDL,NDDL0,IDDL(*)  ,NDOF(*)  ,IKC(*),W_DDL(*),  
     .         ft_n  ,lt_n  ,f_ddl ,l_ddl ,itask
C     REAL
      my_real
     .  lamda,dg(3,*),dgr(3,*),di(3,*),dir(3,*),dd(3,*),ddr(3,*),
     .  l_a,sw2
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .  gi,ui2,uiud,uiut,s
C-----------------------------------------------
         CALL produt_uh(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                 di    ,dir   ,ui2   ,w_ddl  ,f_ddl ,
     .                 l_ddl ,itask )
         CALL produt_uh2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dg    ,dgr   ,di    ,dir   ,uiut  ,
     .                  w_ddl ,f_ddl ,l_ddl ,itask )
         CALL produt_uh2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dd    ,ddr   ,di    ,dir   ,uiud  ,
     .                  w_ddl ,f_ddl ,l_ddl ,itask )
C----------------------
      CALL my_barrier
C---------------------
         gi = sqrt(ui2)
         s = sw2*lamda   
         lamda = ((l_a-gi)*gi-uiud)/(uiut+s)
C
      RETURN
      END
!||====================================================================
!||    b_correct_hp   ../engine/source/implicit/nl_solv.F
!||--- calls      -----------------------------------------------------
!||    imp_smpini     ../engine/source/implicit/imp_solv.F
!||====================================================================
      SUBROUTINE b_correct_hp(NDDL  ,F   ,FEXT  ,DLAMDA )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NDDL
C     REAL
      my_real
     .  dlamda,f(*),fext(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER  F_DDL  ,L_DDL ,ITSK
      INTEGER I
C-----------------------------------------------
!$OMP PARALLEL PRIVATE(ITSK,F_DDL  ,L_DDL,I)
      CALL imp_smpini(itsk   ,f_ddl  ,l_ddl ,nddl )
         DO i =f_ddl  ,l_ddl
                 f(i)=f(i)+dlamda*fext(i)
               ENDDO
!$OMP END PARALLEL
C
      RETURN
      END
!||====================================================================
!||    al_constraint1_hp   ../engine/source/implicit/nl_solv.F
!||--- called by ------------------------------------------------------
!||    nl_solv             ../engine/source/implicit/nl_solv.F
!||--- calls      -----------------------------------------------------
!||    frac_dd_hp          ../engine/source/implicit/integrator.F
!||    integrator2_hp      ../engine/source/implicit/integrator.F
!||    produt_uhp          ../engine/source/implicit/produt_v.F
!||    produt_uhp2         ../engine/source/implicit/produt_v.F
!||====================================================================
      SUBROUTINE al_constraint1_hp(NDDL0 ,NDDL  ,IDDL  ,NDOF  ,IKC   ,
     .                           DD    ,DDR   ,DG    ,DGR   ,DI    ,
     .                           DIR   ,W_DDL ,L_A   ,LAMDA ,SW2   ,
     .                           IER   )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  NDDL,NDDL0,IDDL(*)  ,NDOF(*)  ,IKC(*),W_DDL(*) ,IER,
     .         FT_N  ,LT_N  
C     REAL
      my_real
     .  lamda,dg(3,*),dgr(3,*),di(3,*),dir(3,*),dd(3,*),ddr(3,*),
     .  l_a,sw2
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .  A,B,C,UI2,UIUD,UIUT,UT2,LAM1,LAM2,S,TH1,TH2,FAC
C-----------------------------------------------
         ier=0
         CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                  dg    ,dgr   ,ut2   ,w_ddl  )
C---------------------
         a = ut2+sw2
C-------DD< DI+DD   
         CALL integrator2_hp(iddl,ndof,ikc,dd,ddr,di,dir)
C---------------------
         CALL produt_uhp2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                   dg    ,dgr   ,dd    ,ddr   ,uiut  ,
     .                   w_ddl )
         CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                   dd    ,ddr   ,ui2   ,w_ddl  )
C---------------------
         b= two*(uiut+sw2*lamda)
         c = ui2-l_a*l_a+sw2*lamda
               s= b*b-four*a*c
              IF (s>=0) THEN
                s=sqrt(s)
                lam1=half*(-b+s)/a
                lam2=half*(-b-s)/a
         CALL produt_uhp2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                   dd    ,ddr   ,di    ,dir   ,uiud  ,
     .                   w_ddl )
         CALL produt_uhp2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                   dg    ,dgr   ,di    ,dir   ,uiut  ,
     .                   w_ddl )
C---------------------
         th1= uiud+lam1*uiut
         th2= uiud+lam2*uiut
               IF (th1>zero.AND.th2<zero) THEN
                lamda =lam1
               ELSEIF (th1<zero.AND.th2>zero) THEN
                lamda =lam2
               ELSEIF (th1>zero.AND.th2>zero) THEN
                lamda =min(lam1,lam2)
               ELSE
                ier=1
               ENDIF
        ELSE
               ier=1
              ENDIF
              fac =-one
        CALL frac_dd_hp(iddl,ndof,ikc,dd,ddr,di,dir,fac)
C
      RETURN
      END
!||====================================================================
!||    al_constraint2_hp   ../engine/source/implicit/nl_solv.F
!||--- called by ------------------------------------------------------
!||    nl_solv             ../engine/source/implicit/nl_solv.F
!||--- calls      -----------------------------------------------------
!||    produt_uhp          ../engine/source/implicit/produt_v.F
!||    produt_uhp2         ../engine/source/implicit/produt_v.F
!||====================================================================
      SUBROUTINE al_constraint2_hp(NDDL0 ,NDDL  ,IDDL  ,NDOF  ,IKC   ,
     .                             DD    ,DDR   ,DG    ,DGR   ,DI    ,
     .                             DIR   ,W_DDL ,L_A   ,LAMDA ,SW2   )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  NDDL,NDDL0,IDDL(*)  ,NDOF(*)  ,IKC(*),W_DDL(*)  
C     REAL
      my_real
     .  lamda,dg(3,*),dgr(3,*),di(3,*),dir(3,*),dd(3,*),ddr(3,*),
     .  l_a,sw2
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .  gi,ui2,uiud,uiut,s
C-----------------------------------------------
         CALL produt_uhp(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                   di    ,dir   ,ui2   ,w_ddl  )
         CALL produt_uhp2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                    dg    ,dgr   ,di    ,dir   ,uiut  ,
     .                    w_ddl )
         CALL produt_uhp2(nddl0 ,nddl  ,iddl  ,ndof  ,ikc   ,
     .                    dd    ,ddr   ,di    ,dir   ,uiud  ,
     .                    w_ddl )
C---------------------
         gi = sqrt(ui2)
         s = sw2*lamda   
         lamda = ((l_a-gi)*gi-uiud)/(uiut+s)
C
      RETURN
      END