solve.f

Go to the documentation of this file.

!                   ****************
                    SUBROUTINE solve
!                   ****************
!
     &(x, a,b,tb,cfg,infogr,mesh,aux)
!
!***********************************************************************
! BIEF   V6P2                                   21/08/2010
!***********************************************************************
!
!brief    SOLVES A LINEAR SYSTEM OF THE FORM A X = B
!+
!+            USING ITERATIVE METHODS
!+                WITH POSSIBLE PRECONDITIONING.
!code
!+-----------------------------------------------------------------------
!+                        CHOICE OF THE METHOD
!+-----------------------------------------------------------------------
!+      VALUE          I       MEANING         I
!+-----------------------------------------------------------------------
!+                     I                       I
!+        1            I   CONJUGATE GRADIENT  I
!+                     I                       I
!+        2            I   CONJUGATE RESIDUAL  I
!+                     I                       I
!+        3            I   CONJUGATE GRADIENT  I
!+                     I  ON A NORMAL EQUATION I
!+                     I                       I
!+        4            I     MINIMUM ERROR     I
!+                     I                       I
!+        5            I        SQUARED        I
!+                     I   CONJUGATE GRADIENT  I
!+                     I                       I
!+        6            I        SQUARED        I
!+                     I   CONJUGATE GRADIENT  I
!+                     I        (CGSTAB)       I
!+                     I                       I
!+        7            I       GMRES           I
!+                     I                       I
!+        8            I     DIRECT : YSMP     I
!+                     I                       I
!+        9            I     DIRECT : MUMPS    I
!+                     I                       I
!+-----------------------------------------------------------------------
!
!warning  FOR SOME PRECONDITIONING ALGORITHMS, MATRIX A
!+            CAN BE MODIFIED
!code
!+                        PRECONDITIONING
!+
!+     NOTES     : SOME PRECONDITIONING ALGORITHMS CAN BE ADDED
!+                 (DIAGONAL 2, 3 OR 5 WITH THE OTHERS).
!+                 THAT'S WHY PRIME NUMBERS WERE USED TO CHARACTERISE
!+                 THE PRECONDITIONING ALGORITHMS.
!+                 SHOULD THE USER WANT TO CUMULATE THE EFFECTS, HE/SHE
!+                 SHOULD GIVE THE PRODUCT OF THE CORRESPONDING PRIMES
!+
!+-----------------------------------------------------------------------
!+        VALUE        I                  MEANING
!+-----------------------------------------------------------------------
!+        0 OR 1       I  NO PRECONDITIONING
!+                     I
!+        2            I  DIAGONAL PRECONDITIONING WITH THE MATRIX
!+                     I  DIAGONAL
!+                     I
!+        3            I  BLOCK DIAGONAL PRECONDITIONING
!+                     I
!+        5            I  DIAGONAL PRECONDITIONING WITH THE ABSOLUTE
!+                     I  VALUE OF THE MATRIX DIAGONAL
!+                     I
!+        7            I  CROUT'S PRECONDITIONING BY ELEMENT
!+                     I
!+       11            I  GAUSS-SEIDEL'S PRECONDITIONING BY ELEMENT
!+                     I
!+       13            I  PRECONDITIONING MATRIX SUPPLIED BY THE USER
!+                     I
!+-----------------------------------------------------------------------
!
!history  J-M HERVOUET (LNHE)
!+        18/02/08
!+        V5P9
!+
!
!history  C. DENIS (SINETICS)
!+        19/03/10
!+        V6P0
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!history  J-M HERVOUET (LNHE)
!+        07/12/2011
!+        V6P2
!+   Call to preverseg and preverebe modified.
!+
!
!history  J.PARISI (HRW)
!+        9/08/2012
!+        V6P2
!+   Call to SD_SOLVE_4 modified.
!history  R.NHEILI (Univerte de Perpignan, DALI)
!+        24/02/2016
!+        V7P3
!+        COMPENSATED INTERFACE NODE ASSEMBLY (MODASS=3)
!
!history  J-M HERVOUET (EDF LAB, LNHE)
!+        23/02/2015
!+        V7P1
!+   Arguments added to PRE4_MUMPS, to correct a bug when it calls
!+   SD_FABSG4.
!
!history  J,RIEHME (ADJOINTWARE)
!+        November 2016
!+        V7P2
!+   Replaced EXTERNAL statements to parallel functions / subroutines
!+   by the INTERFACE_PARALLEL
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| A              |-->| MATRIX OF THE SYSTEM (OR BLOCK OF MATRICES)
!| AUX            |-->| MATRIX FOR PRECONDITIONING.
!| B              |-->| RIGHT-HAND SIDE OF THE SYSTEM
!| CFG            |-->| STRUCTURE OF SOLVER CONFIGURATION
!|                |   | CFG%KRYLOV IS USED ONLY IF CFG%SLV = 7 (GMRES)
!| INFOGR         |-->| IF YES, PRINT A LOG.
!| MESH           |-->| MESH STRUCTURE.
!| TB             |-->| BLOCK OF VECTORS WITh AT LEAST
!|                |   | MAX(7,2+2*CFG%KRYLOV)*S VECTORS, S IS 1
!|                |   | IF A IS A MATRIX, 2 IF A BLOCK OF 4 MATRICES
!|                |   | AND 3 IF A BLOCK OF 9.
!| X              |<->| INITIAL VALUE, THEN SOLUTION
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_solve => solve
      USE declarations_telemac, ONLY : modass
      USE declarations_telemac, ONLY : tbb, bb, bx, first_solve
      USE declarations_special
!
      USE interface_parallel, ONLY : p_max
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      TYPE(slvcfg), INTENT(INOUT) :: CFG
!
!     STRUCTURES OF VECTORS OR BLOCKS OF VECTORS
!
      TYPE(bief_obj), TARGET, INTENT(INOUT) :: X,B
      TYPE(bief_obj), INTENT(INOUT)         :: TB
!
!     STRUCTURES OF MATRIX OR BLOCK OF MATRICES
!
      TYPE(bief_obj), INTENT(INOUT) :: A,AUX
!
      LOGICAL, INTENT(IN) :: INFOGR
!
!     MESH STRUCTURE
!
      TYPE(bief_mesh), INTENT(INOUT) :: MESH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER PRESTO,IG,LV,S,NBL,I
      INTEGER IT1,IT2,IT3,IT4,IT5,IT6,IT7,IBL1,IBL2,K,IAD,ITB,ITBB
!
      LOGICAL DIADON,PREXSM
      INTEGER NPOIN_TOT
!
      INTRINSIC max
!
!-----------------------------------------------------------------------
!
!     STRUCTURES OF BLOCKS OF WORKING ARRAYS
!
      TYPE(bief_obj), POINTER :: PB,PX
!
!-----------------------------------------------------------------------
!
!  ALLOCATES THE BLOCK OF BLOCKS TBB AND THE BLOCKS IN TBB
!
      IF(first_solve) THEN
        CALL allblo(tbb,'TBB   ')
        CALL allblo(bb ,'BB    ')
        CALL allblo(bx ,'BX    ')
        first_solve=.false.
      ENDIF
      nbl = 7
      IF(cfg%SLV.EQ.7) nbl = max(nbl,4+2*cfg%KRYLOV)
      IF(nbl.GT.tbb%N) THEN
        IF(tbb%N.NE.0) THEN
          DO i=1,tbb%N
            DEALLOCATE(tbb%ADR(i)%P%ADR)
            NULLIFY(tbb%ADR(i)%P%ADR)
            DEALLOCATE(tbb%ADR(i)%P)
            NULLIFY(tbb%ADR(i)%P)
          ENDDO
        ENDIF
        tbb%N=0
        CALL allblo_in_block(tbb,nbl,'BL    ')
      ENDIF
!
!-----------------------------------------------------------------------
!
      lv = mesh%LV
!
!-----------------------------------------------------------------------
!
!  VARIOUS TYPES OF SOLVED SYSTEMS S = 0 : NORMAL MATRIX
!                                  S = 1 : 1 MATRICE IN A BLOCK
!                                  S = 2 : 4 MATRICES IN A BLOCK
!                                  S = 3 : 9 MATRICES IN A BLOCK
!
      IF(a%TYPE.EQ.3) THEN
        s = 0
        bb%N = 0
        bx%N = 0
        CALL addblo(bb,b)
        CALL addblo(bx,x)
        px => bx
        pb => bb
      ELSEIF(a%TYPE.EQ.4) THEN
        IF(a%N.EQ.1) THEN
          s = 1
        ELSEIF(a%N.EQ.4) THEN
          s = 2
        ELSEIF(a%N.EQ.9) THEN
          s = 3
        ENDIF
        px => x
        pb => b
      ENDIF
!
!----------------
!  DIRECT SOLVERS
!  ---> YSMP
!----------------
!
      IF(cfg%SLV.EQ.8) THEN
!
        IF(ncsize.GT.1) THEN
          WRITE(lu,2019)
2019      FORMAT(1x,'USE THE PARALLEL DIRECT SOLVER MUMPS,',/,1x,
     &              'SOLVER = 9',///)
          CALL plante(1)
          stop
        ENDIF
!
        IF(s.EQ.0) THEN
          IF(a%TYPEXT.NE.'S'.AND.a%TYPEXT.NE.'Q') THEN
            WRITE(lu,*) 'SOLVE (BIEF): OFF-DIAGONAL TERMS'
            WRITE(lu,*) '              OF TYPE ',a%TYPEXT
            WRITE(lu,*) '              NOT IMPLEMENTED'
            CALL plante(1)
            stop
          ENDIF
          CALL sd_solve_1(a%D%DIM1,mesh%NSEG,mesh%GLOSEG%I,
     &                    mesh%GLOSEG%DIM1,
     &                    a%D%R,a%X%R,x%R,b%R,infogr,a%TYPEXT)
        ELSEIF(s.EQ.1) THEN
          IF(a%ADR(1)%P%TYPEXT.NE.'S'.AND.a%ADR(1)%P%TYPEXT.NE.'Q') THEN
            WRITE(lu,*) 'SOLVE (BIEF): DIRECT SOLVER FOR SYMMETRIC'
            WRITE(lu,*) '              SYSTEMS ONLY'
            CALL plante(1)
            stop
          ENDIF
          CALL sd_solve_1(a%ADR(1)%P%D%DIM1,mesh%NSEG,mesh%GLOSEG%I,
     &                    mesh%GLOSEG%DIM1,
     &                    a%ADR(1)%P%D%R,a%ADR(1)%P%X%R,x%ADR(1)%P%R,
     &                    b%ADR(1)%P%R,infogr,a%ADR(1)%P%TYPEXT)
        ELSEIF(s.EQ.2) THEN
          CALL sd_solve_4(mesh%NPOIN,mesh%NSEG,mesh%GLOSEG%I,
     &                    a%ADR(1)%P%D%R,a%ADR(2)%P%D%R,
     &                    a%ADR(3)%P%D%R,a%ADR(4)%P%D%R,
     &                    a%ADR(1)%P%X%R,a%ADR(2)%P%X%R,
     &                    a%ADR(3)%P%X%R,a%ADR(4)%P%X%R,
     &                    x%ADR(1)%P%R,x%ADR(2)%P%R,
     &                    b%ADR(1)%P%R,b%ADR(2)%P%R,infogr,
     &                    a%ADR(1)%P%TYPEXT,a%ADR(2)%P%TYPEXT,
     &                    a%ADR(3)%P%TYPEXT,a%ADR(4)%P%TYPEXT)
!       ELSEIF(S.EQ.3) THEN
        ELSE
          WRITE(lu,401) s
401       FORMAT(1x,'SOLVE (BIEF): S=',1i6,' CASE NOT IMPLEMENTED')
          CALL plante(1)
          stop
        ENDIF
        RETURN
!
      ELSEIF(cfg%SLV.EQ.9) THEN
!
!----------------
!  DIRECT SOLVERS
!  ---> MUMPS
!----------------
!
        IF(ncsize.LT.1) THEN
          WRITE(lu,3019)
3019      FORMAT(1x,'MUMPS ARE NOT AVAILABLE FOR SEQUENTIAL RUNS,',/,1x,
     &         'USE SEQUENITAL DIRECT SOLVER (SOLVER = 8) ',///)
          CALL plante(1)
          stop
        ENDIF
!
!       COMPUTING THE NUMBER OF POINTS IN THE MESH BEFORE PARTITIONING
!
        IF(ncsize.GT.1) THEN
          npoin_tot=0
          DO i=1,mesh%NPOIN
            npoin_tot=max(mesh%KNOLG%I(i),npoin_tot)
          ENDDO
          npoin_tot=p_max(npoin_tot)
        ELSE
          npoin_tot=mesh%NPOIN
        ENDIF
!
        IF(s.EQ.0) THEN
          WRITE(lu,402) s
402       FORMAT(1x,'SOLVE (BIEF): S=',1i6,1x,
     &              'CASE NOT YET MPLEMENTED FOR MUMPS')
          CALL plante(1)
          stop
        ELSEIF(s.EQ.1) THEN
          WRITE(lu,4011) s
 4011     FORMAT(1x,'SOLVE (BIEF): S=',1i6,' CASE NOT YET
     &           IMPLEMENTED FOR MUMPS')
          CALL plante(1)
          stop
        ELSEIF(s.EQ.2) THEN
          CALL pre4_mumps(mesh%NPOIN,mesh%NSEG,mesh%GLOSEG%I,
     &                    a%ADR(1)%P%D%R,a%ADR(2)%P%D%R,
     &                    a%ADR(3)%P%D%R,a%ADR(4)%P%D%R,
     &                    a%ADR(1)%P%X%R,a%ADR(2)%P%X%R,
     &                    a%ADR(3)%P%X%R,a%ADR(4)%P%X%R,
     &                    x%ADR(1)%P%R,x%ADR(2)%P%R,
     &                    b%ADR(1)%P%R,b%ADR(2)%P%R,
     &                    a%ADR(1)%P%TYPEXT,a%ADR(2)%P%TYPEXT,
     &                    a%ADR(3)%P%TYPEXT,a%ADR(4)%P%TYPEXT,
     &                    mesh%KNOLG%I,npoin_tot)
        ELSE
          WRITE(lu,401) s
          CALL plante(1)
          stop
        ENDIF
        RETURN
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      presto = cfg%PRECON
      IF(cfg%PRECON.EQ.0) cfg%PRECON = 1
!
!-----------------------------------------------------------------------
!
!  MANAGES WORKING ARRAYS : ITB --> NEXT AVAILABLE VECTOR
!
!  ITB  --> NEXT AVAILABLE VECTOR IN TB
!  ITBB --> NEXT AVAIALBLE BLOCK  IN TBB
!
      itb  = 1
      itbb = 1
!
!  ALLOCATES TWO WORKING BLOCKS WITH, EITHER A VECTOR,
!  OR A BLOCK OF VECTORS (CASE WHERE S IS 0).
!  THESE TWO BLOCKS ARE COMMON TO ALL THE METHODS.
!
!     FOR THE PRECONDITIONING MATRICES
      IF(3*(cfg%PRECON/3).EQ.cfg%PRECON) THEN
!       BLOCK DIAGONAL PRECONDITIONING : S**2 DIAGONALS
        CALL solaux(it1, tb,tbb,itb,itbb,s**2)
      ELSE
!       OTHER : S DIAGONALS
        CALL solaux(it1, tb,tbb,itb,itbb,s)
      ENDIF
!
      CALL solaux(it2, tb,tbb,itb,itbb,s)
!
      IF(cfg%SLV.EQ.7) THEN
!       SPECIAL GMRES : ARRAYS DEPENDING ON THE SIZE OF KRYLOV
!                       TBB(IBL1) : BLOCK OF CFG%KRYLOV VECTORS
!                       OR BLOCK OF CFG%KRYLOV BLOCKS OF S VECTORS
!                       TBB(IBL2) : IDEM
!
        ibl1=itbb
        itbb = itbb + 1
        ibl2=itbb
        itbb = itbb + 1
        DO i=1,tbb%ADR(ibl1)%P%N
          NULLIFY(tbb%ADR(ibl1)%P%ADR(i)%P)
        ENDDO
        tbb%ADR(ibl1)%P%N=0
        DO i=1,tbb%ADR(ibl2)%P%N
          NULLIFY(tbb%ADR(ibl2)%P%ADR(i)%P)
        ENDDO
        tbb%ADR(ibl2)%P%N=0
        DO k=1,cfg%KRYLOV
          CALL solaux(iad, tb,tbb,itb,itbb,s)
          CALL addblo(tbb%ADR(ibl1)%P,tbb%ADR(iad)%P)
          CALL solaux(iad, tb,tbb,itb,itbb,s)
          CALL addblo(tbb%ADR(ibl2)%P,tbb%ADR(iad)%P)
        ENDDO ! K
!       AVOIDS A WARNING FROM THE INTEL COMPILER
        it3=-1
        it4=-1
        it5=-1
        it6=-1
        it7=-1
      ELSE
!       OTHER METHODS (COULD SOMETIMES NOT ALLOCATE IT6 OR IT7)
        CALL solaux(it3, tb,tbb,itb,itbb,s)
        CALL solaux(it4, tb,tbb,itb,itbb,s)
        CALL solaux(it5, tb,tbb,itb,itbb,s)
        CALL solaux(it6, tb,tbb,itb,itbb,s)
        CALL solaux(it7, tb,tbb,itb,itbb,s)
!       AVOIDS A WARNING FROM THE CRAY COMPILER
        ibl1=1
        ibl2=1
!
      ENDIF
!
!     CROUT'S PRECONDITIONING : REQUIRES A PRELIMINARY PRECONDITIONING
!     THAT SETS DIAGONALS TO 1.
!     THE GRADIENT WILL BE DISTINGUISHED FROM THE RESIDUE (POINTER IG)
!
      IF(  7*(cfg%PRECON/ 7).EQ.cfg%PRECON.OR.
     &    11*(cfg%PRECON/11).EQ.cfg%PRECON.OR.
     &    13*(cfg%PRECON/13).EQ.cfg%PRECON.OR.
     &    17*(cfg%PRECON/17).EQ.cfg%PRECON     ) THEN
        ig=it6
        IF(2*(cfg%PRECON/2).NE.cfg%PRECON.AND.
     &     3*(cfg%PRECON/3).NE.cfg%PRECON) THEN
!         SELECTS DIAGONAL
          cfg%PRECON=2*cfg%PRECON
        ENDIF
      ELSE
!       NOTE IT5 =-1 IF CFG%SLV.EQ.7 BUT IN THIS CASE IG IS USELESS
        ig=it5
      ENDIF
!
!  END OF: MANAGES THE WORKING ARRAYS
!
!-----------------------------------------------------------------------
!                               -1/2      -1/2  1/2        -1/2
!  DIAGONAL PRECONDITIONINGS : D     A  D      D     X  = D      B
!
      diadon = .false.
      prexsm = .true.
!
      IF(3*(cfg%PRECON/3).EQ.cfg%PRECON.AND.(s.EQ.2.OR.s.EQ.3)) THEN
!       BLOCK-DIAGONAL PRECONDITIONING (4 OR 9 MATRICES)
        CALL prebdt(x,a,b,tbb%ADR(it1)%P,mesh,prexsm,diadon,s)
!       DOES NOT MOFIFY D11, D22, D33 WHEN PRECDT IS CALLED
        diadon = .true.
      ENDIF
!
      IF(2*(cfg%PRECON/2).EQ.cfg%PRECON.OR.
     &   3*(cfg%PRECON/3).EQ.cfg%PRECON.OR.
     &   5*(cfg%PRECON/5).EQ.cfg%PRECON) THEN
        CALL precdt(x,a,b,tbb%ADR(it1)%P,mesh,
     &              cfg%PRECON,prexsm,diadon,s)
      ENDIF
!
!-----------------------------------------------------------------------
!
!  BUILDS THE PRECONDITIONING MATRICES:
!
      IF(7*(cfg%PRECON/7).EQ.cfg%PRECON) THEN
        CALL dcpldu(aux,a,mesh,.true.,lv)
      ELSEIF(11*(cfg%PRECON/11).EQ.cfg%PRECON) THEN
        CALL gsebe(aux,a,mesh)
      ELSEIF(13*(cfg%PRECON/13).EQ.cfg%PRECON) THEN
!       DOES NOTHING, AUX IS SUPPLIED BY THE SUBROUTINE CALLING
      ELSEIF(17*(cfg%PRECON/17).EQ.cfg%PRECON) THEN
        IF(cfg%SLV.NE.1.AND.cfg%SLV.NE.2) THEN
          WRITE(lu,*) 'PRECONDITIONING 17'
          WRITE(lu,*) 'NOT IMPLEMENTED FOR SOLVER ',cfg%SLV
          CALL plante(1)
          stop
        ENDIF
        IF(aux%TYPE.NE.3) THEN
          WRITE(lu,*) 'PRECONDITIONING 17'
          WRITE(lu,*) 'NOT IMPLEMENTED FOR BLOCKS OF MATRICES'
          CALL plante(1)
          stop
        ENDIF
!
        IF(aux%STO.EQ.1) THEN
          CALL preverebe(aux%X%R,a%X%R,a%TYPDIA,a%TYPEXT,
     &                   mesh%IKLE%I,mesh%NPOIN,mesh%NELEM,
     &                   mesh%NELMAX,mesh,mesh%TYPELM)
        ELSE
          CALL preverseg(aux%X%R,a%D%R,a%X%R,a%TYPDIA,a%TYPEXT,
     &                   mesh%NPOIN,mesh,mesh%NSEG,mesh%TYPELM)
        ENDIF
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!  PARALLEL MODE: SECOND MEMBER
!
      IF(ncsize.GT.1) THEN
        IF (modass .EQ. 1) THEN
          CALL parcom(b,2,mesh)
        ELSEIF (modass .EQ. 3) THEN
          CALL parcom_comp(b,b%E,2,mesh)
        ENDIF
      ENDIF
      IF (modass .EQ.3)THEN
        b%R=b%R+b%E
        b%E=0.d0
      ENDIF
!
!-----------------------------------------------------------------------
!
!  SOLVES THE LINEAR SYSTEM:
!
      IF(cfg%SLV.EQ.1) THEN
!
!       CONJUGATE GRADIENT
!
        CALL gracjg(px, a,pb, mesh,
     &              tbb%ADR(it2)%P,tbb%ADR(it3)%P,
     &              tbb%ADR(it5)%P,tbb%ADR(ig)%P,
     &              cfg,infogr,aux)
!
      ELSEIF(cfg%SLV.EQ.2) THEN
!
!       CONJUGATE RESIDUAL
!
        CALL rescjg(px, a,pb, mesh,
     &              tbb%ADR(it2)%P,tbb%ADR(it3)%P,
     &              tbb%ADR(it4)%P,tbb%ADR(it5)%P,
     &              tbb%ADR(ig)%P,
     &              cfg,infogr,aux)
!
      ELSEIF(cfg%SLV.EQ.3) THEN
!
!       NORMAL EQUATION
!
        CALL equnor(px, a,pb, mesh,
     &              tbb%ADR(it2)%P,tbb%ADR(it3)%P,
     &              tbb%ADR(it4)%P,tbb%ADR(it5)%P,
     &              tbb%ADR(ig)%P,
     &              cfg,infogr,aux)
!
      ELSEIF(cfg%SLV.EQ.4) THEN
!
!       MINIMUM ERROR
!
        CALL errmin(px, a,pb, mesh,
     &              tbb%ADR(it2)%P,tbb%ADR(it3)%P,tbb%ADR(it5)%P,
     &              tbb%ADR(ig)%P,
     &              cfg,infogr,aux)
!
      ELSEIF(cfg%SLV.EQ.5) THEN
!
!       SQUARED CONJUGATE GRADIENT
!
        CALL cgsqua(px, a,pb, mesh,
     &              tbb%ADR(it2)%P,tbb%ADR(it3)%P,tbb%ADR(it4)%P,
     &              tbb%ADR(it5)%P,tbb%ADR(it6)%P,tbb%ADR(it7)%P,
     &              cfg,infogr)
!
      ELSEIF(cfg%SLV.EQ.6) THEN
!
!       CGSTAB
!
        CALL cgstab(px, a,pb, mesh,
     &              tbb%ADR(it2)%P,tbb%ADR(it3)%P,tbb%ADR(it4)%P,
     &              tbb%ADR(it5)%P,tbb%ADR(it6)%P,tbb%ADR(it7)%P,
     &              cfg,infogr,aux)
!
      ELSEIF(cfg%SLV.EQ.7) THEN
!
!       GENERALISED MINIMUM RESIDUAL
!
        CALL gmres(px, a,pb,mesh,
     &             tbb%ADR(it2)%P,tbb%ADR(ibl1)%P,tbb%ADR(ibl2)%P,
     &             cfg,infogr,aux)
!
!
      ELSE
!
        WRITE(lu,400) cfg%SLV
400     FORMAT(1x,'SOLVE (BIEF) :',1i6,' METHOD NOT AVAILABLE :')
        CALL plante(1)
        stop
!
      ENDIF
!
!-----------------------------------------------------------------------
!
!  INVERSES THE CHANGE IN VARIABLE IF PRECONDITIONING
!                                                DIAGONAL
!                                                DIAGONAL-BLOCK
!
      IF(2*(cfg%PRECON/2).EQ.cfg%PRECON.OR.
     &   3*(cfg%PRECON/3).EQ.cfg%PRECON.OR.
     &   5*(cfg%PRECON/5).EQ.cfg%PRECON    ) THEN
        CALL os('X=XY    ', x=px, y=tbb%ADR(it1)%P)
      ENDIF
!
      IF(3*(cfg%PRECON/3).EQ.cfg%PRECON.AND.(s.EQ.2.OR.s.EQ.3)) THEN
        CALL um1x(x,tbb%ADR(it1)%P,s)
      ENDIF
!
!-----------------------------------------------------------------------
!
      cfg%PRECON = presto
!
!-----------------------------------------------------------------------
!
      RETURN
      END