sfac_scalings.F

Go to the documentation of this file.

C
C  This file is part of MUMPS 5.5.1, released
C  on Tue Jul 12 13:17:24 UTC 2022
C
C
C  Copyright 1991-2022 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria,
C  Mumps Technologies, University of Bordeaux.
C
C  This version of MUMPS is provided to you free of charge. It is
C  released under the CeCILL-C license 
C  (see doc/CeCILL-C_V1-en.txt, doc/CeCILL-C_V1-fr.txt, and
C  https://cecill.info/licences/Licence_CeCILL-C_V1-en.html)
C
       SUBROUTINE smumps_fac_a(N, NZ8, NSCA, 
     &      ASPK, IRN, ICN, COLSCA, ROWSCA, WK, LWK8, WK_REAL,
     &      LWK_REAL, ICNTL, INFO)
       IMPLICIT NONE
      INTEGER N, NSCA
      INTEGER(8), INTENT(IN) :: NZ8
      INTEGER IRN(NZ8), ICN(NZ8)
      INTEGER ICNTL(60), INFO(80)
      REAL, INTENT(IN) :: ASPK(NZ8)
      REAL COLSCA(*), ROWSCA(*)
      INTEGER(8), INTENT(IN) :: LWK8
      INTEGER LWK_REAL
      REAL    WK(LWK8)
      REAL WK_REAL(LWK_REAL)
      INTEGER MPG,LP
      INTEGER IWNOR
      INTEGER I
      LOGICAL PROK
      REAL ONE
      parameter( one = 1.0e0 )
      lp      = icntl(1)
      mpg     = icntl(2)
      mpg    = icntl(3)
      prok   = ((mpg.GT.0).AND.(icntl(4).GE.2))
      IF (prok) THEN
        WRITE(mpg,101)
      ELSE
        mpg = 0
      ENDIF
 101    FORMAT(/' ****** SCALING OF ORIGINAL MATRIX '/)
        IF (nsca.EQ.1) THEN
         IF (prok)
     &    WRITE (mpg,*) ' DIAGONAL SCALING '
        ELSEIF (nsca.EQ.3) THEN
         IF (prok)
     &   WRITE (mpg,*) ' COLUMN SCALING'
        ELSEIF (nsca.EQ.4) THEN
         IF (prok)
     &   WRITE (mpg,*) ' ROW AND COLUMN SCALING (1 Pass)'
       ENDIF
        DO 10 i=1,n
            colsca(i) = one
            rowsca(i) = one
 10     CONTINUE
        IF (5*n.GT.lwk_real) GOTO 410
        iwnor = 1
          IF (nsca.EQ.1) THEN
            CALL smumps_fac_v(n,nz8,aspk,irn,icn,
     &        colsca,rowsca,mpg)
          ELSEIF (nsca.EQ.3) THEN
            CALL smumps_fac_y(n,nz8,aspk,irn,icn,wk_real(iwnor),
     &      colsca, mpg)
          ELSEIF (nsca.EQ.4) THEN
            CALL smumps_rowcol(n,nz8,irn,icn,aspk,
     &      wk_real(iwnor),wk_real(iwnor+n),colsca,rowsca,mpg)
          ENDIF
      GOTO 500
 410  info(1) = -5
      info(2) = 5*n-lwk_real
      IF ((lp.GT.0).AND.(icntl(4).GE.1))
     & WRITE(lp,*) '*** ERROR: Not enough space to scale matrix'
      GOTO 500
 500  CONTINUE
      RETURN
      END SUBROUTINE smumps_fac_a
      SUBROUTINE smumps_rowcol(N,NZ8,IRN,ICN,VAL,
     &    RNOR,CNOR,COLSCA,ROWSCA,MPRINT)
      INTEGER,    INTENT(IN) :: N
      INTEGER(8), INTENT(IN) :: NZ8
      REAL    VAL(NZ8)
      REAL    RNOR(N),CNOR(N)
      REAL    COLSCA(N),ROWSCA(N)
      REAL    CMIN,CMAX,RMIN,ARNOR,ACNOR
      INTEGER IRN(NZ8), ICN(NZ8)
      REAL    VDIAG
      INTEGER MPRINT
      INTEGER I,J
      INTEGER(8) :: K8
      REAL ZERO, ONE
      parameter(zero=0.0e0, one=1.0e0)
      DO 50 j=1,n
       cnor(j)   = zero
       rnor(j)   = zero
  50  CONTINUE
      DO 100 k8=1_8,nz8
          i = irn(k8)
          j = icn(k8)
          IF ((i.LE.0).OR.(i.GT.n).OR.
     &        (j.LE.0).OR.(j.GT.n)) GOTO 100
            vdiag = abs(val(k8))
            IF (vdiag.GT.cnor(j)) THEN
              cnor(j) =     vdiag
            ENDIF
            IF (vdiag.GT.rnor(i)) THEN
              rnor(i) =     vdiag
            ENDIF
 100   CONTINUE
      IF (mprint.GT.0) THEN
       cmin = cnor(1)
       cmax = cnor(1)
       rmin = rnor(1)
       DO 111 i=1,n
        arnor = rnor(i)
        acnor = cnor(i)
        IF (acnor.GT.cmax) cmax=acnor
        IF (acnor.LT.cmin) cmin=acnor
        IF (arnor.LT.rmin) rmin=arnor
 111   CONTINUE
       WRITE(mprint,*) '**** STAT. OF MATRIX PRIOR ROW&COL SCALING'
       WRITE(mprint,*) ' MAXIMUM NORM-MAX OF COLUMNS:',cmax
       WRITE(mprint,*) ' MINIMUM NORM-MAX OF COLUMNS:',cmin
       WRITE(mprint,*) ' MINIMUM NORM-MAX OF ROWS   :',rmin
      ENDIF
      DO 120 j=1,n
       IF (cnor(j).LE.zero) THEN
         cnor(j)   = one
       ELSE
         cnor(j)   = one / cnor(j)
       ENDIF
 120  CONTINUE
      DO 130 j=1,n
       IF (rnor(j).LE.zero) THEN
         rnor(j)   = one
       ELSE
         rnor(j)   = one / rnor(j)
       ENDIF
 130  CONTINUE
       DO 110 i=1,n
        rowsca(i) = rowsca(i) * rnor(i)
        colsca(i) = colsca(i) * cnor(i)
 110   CONTINUE
      IF (mprint.GT.0)
     &  WRITE(mprint,*) ' END OF SCALING BY MAX IN ROW AND COL'
      RETURN
      END SUBROUTINE smumps_rowcol
      SUBROUTINE smumps_fac_y(N,NZ8,VAL,IRN,ICN,
     &       CNOR,COLSCA,MPRINT)
      INTEGER,    INTENT(IN)  :: N
      INTEGER(8), INTENT(IN)  :: NZ8
      REAL,    INTENT(IN)  :: VAL(NZ8)
      REAL,       INTENT(OUT) :: CNOR(N)
      REAL,       INTENT(INOUT) :: COLSCA(N)
      INTEGER,    INTENT(IN)  :: IRN(NZ8), ICN(NZ8)
      INTEGER,    INTENT(IN)  :: MPRINT
      REAL VDIAG
      INTEGER I,J
      INTEGER(8) :: K8
      REAL ZERO, ONE
      parameter(zero=0.0e0,one=1.0e0)
      DO 10 j=1,n
       cnor(j)   = zero
  10  CONTINUE
      DO 100 k8=1_8,nz8
        i = irn(k8)
        j = icn(k8)
        IF ((i.LE.0).OR.(i.GT.n).OR.
     &      (j.LE.0).OR.(j.GT.n)) GOTO 100
        vdiag = abs(val(k8))
        IF (vdiag.GT.cnor(j)) THEN
           cnor(j) =     vdiag
        ENDIF
 100  CONTINUE
      DO 110 j=1,n
       IF (cnor(j).LE.zero) THEN
         cnor(j)   = one
       ELSE
         cnor(j)   = one/cnor(j)
       ENDIF
 110  CONTINUE
       DO 215 i=1,n
        colsca(i) = colsca(i) * cnor(i)
 215   CONTINUE
      IF (mprint.GT.0) WRITE(mprint,*) ' END OF COLUMN SCALING'
      RETURN
      END SUBROUTINE smumps_fac_y
      SUBROUTINE smumps_fac_v(N,NZ8,VAL,IRN,ICN,
     &      COLSCA,ROWSCA,MPRINT)
      INTEGER   , INTENT(IN)  :: N
      INTEGER(8), INTENT(IN)  :: NZ8
      REAL   , INTENT(IN)  :: VAL(NZ8)
      REAL      , INTENT(OUT) :: ROWSCA(N),COLSCA(N)
      INTEGER   , INTENT(IN)  :: IRN(NZ8),ICN(NZ8)
      INTEGER   , INTENT(IN)  :: MPRINT
      REAL       :: VDIAG
      INTEGER    :: I,J
      INTEGER(8) :: K8
      INTRINSIC sqrt
      REAL ZERO, ONE
      parameter(zero=0.0e0, one=1.0e0)
      DO 10 i=1,n
       rowsca(i)   = one
  10  CONTINUE
      DO 100 k8=1_8,nz8
          i = irn(k8)
          IF ((i.GT.n).OR.(i.LE.0)) GOTO 100
          j = icn(k8)
          IF (i.EQ.j) THEN
            vdiag = abs(val(k8))
            IF (vdiag.GT.zero) THEN
              rowsca(j) = one/(sqrt(vdiag))
            ENDIF
          ENDIF
 100   CONTINUE
       DO 110 i=1,n
        colsca(i) = rowsca(i)
 110   CONTINUE
      IF (mprint.GT.0) WRITE(mprint,*) ' END OF DIAGONAL SCALING'
      RETURN
      END SUBROUTINE smumps_fac_v
      SUBROUTINE smumps_fac_x(NSCA,N,NZ8,IRN,ICN,VAL,
     &    RNOR,ROWSCA,MPRINT)
      INTEGER,    INTENT(IN) :: N, NSCA
      INTEGER(8), INTENT(IN) :: NZ8
      INTEGER,    INTENT(IN) :: IRN(NZ8), ICN(NZ8)
      REAL VAL(NZ8)
      REAL RNOR(N)
      REAL ROWSCA(N)
      INTEGER MPRINT
      REAL VDIAG
      INTEGER I,J
      INTEGER(8) :: K8
      REAL, PARAMETER :: ZERO = 0.0e0
      REAL, PARAMETER :: ONE  = 1.0e0
      DO 50 j=1,n
       rnor(j)   = zero
  50  CONTINUE
      DO 100 k8=1_8,nz8
          i = irn(k8)
          j = icn(k8)
          IF ((i.LE.0).OR.(i.GT.n).OR.
     &        (j.LE.0).OR.(j.GT.n)) GOTO 100
            vdiag = abs(val(k8))
            IF (vdiag.GT.rnor(i)) THEN
              rnor(i) =  vdiag
            ENDIF
 100   CONTINUE
      DO 130 j=1,n
       IF (rnor(j).LE.zero) THEN
         rnor(j)   = one
       ELSE
         rnor(j)   = one/rnor(j)
       ENDIF
 130  CONTINUE
      DO 110 i=1,n
        rowsca(i) = rowsca(i)* rnor(i)
 110  CONTINUE
      IF ( (nsca.EQ.4) .OR. (nsca.EQ.6) ) THEN
        DO 150 k8 = 1_8, nz8
          i   = irn(k8)
          j   = icn(k8)
          IF (min(i,j).LT.1 .OR. i.GT.n .OR. j.GT.n) GOTO 150
          val(k8) = val(k8) * rnor(i)
 150    CONTINUE
      ENDIF
      IF (mprint.GT.0)
     &  WRITE(mprint,'(A)') '  END OF ROW SCALING'
      RETURN
      END SUBROUTINE smumps_fac_x
      SUBROUTINE smumps_anorminf( id,  ANORMINF, LSCAL, 
     &            EFF_SIZE_SCHUR )
      USE smumps_struc_def
      IMPLICIT NONE
      include 'mpif.h'
      INTEGER MASTER, IERR
      parameter( master = 0 )
      TYPE(smumps_struc), TARGET :: id
      REAL, INTENT(OUT) :: ANORMINF
      LOGICAL, INTENT(IN)  :: LSCAL
      INTEGER, INTENT(IN)  :: EFF_SIZE_SCHUR
      INTEGER, DIMENSION (:), POINTER :: KEEP,INFO
      INTEGER(8), DIMENSION (:), POINTER :: KEEP8
      LOGICAL :: I_AM_SLAVE
      REAL DUMMY(1)
      REAL ZERO
      parameter( zero = 0.0e0)
      REAL, ALLOCATABLE :: SUMR(:), SUMR_LOC(:)
      INTEGER :: allocok, MTYPE, I
      info =>id%INFO
      keep =>id%KEEP
      keep8 =>id%KEEP8
      i_am_slave = ( id%MYID .ne. master  .OR.
     &             ( id%MYID .eq. master .AND.
     &               keep(46) .eq. 1 ) )
      IF (id%MYID .EQ. master) THEN
       ALLOCATE( sumr( id%N ), stat =allocok )
       IF (allocok .GT.0 ) THEN
        id%INFO(1)=-13
        id%INFO(2)=id%N
        RETURN
       ENDIF
      ENDIF
      IF ( keep(54) .eq. 0 ) THEN
          IF (id%MYID .EQ. master) THEN
            IF (keep(55).EQ.0) THEN
             IF (.NOT.lscal) THEN
              CALL smumps_sol_x(id%A(1),
     &          id%KEEP8(28), id%N,
     &          id%IRN(1), id%JCN(1),
     &          sumr, keep(1),keep8(1),
     &          eff_size_schur, id%SYM_PERM(1) )
             ELSE
              CALL smumps_scal_x(id%A(1),
     &          id%KEEP8(28), id%N,
     &          id%IRN(1), id%JCN(1), 
     &          sumr, keep(1), keep8(1),
     &          id%COLSCA(1), 
     &          eff_size_schur, id%SYM_PERM(1) )
             ENDIF
            ELSE
             mtype = 1
             IF (.NOT.lscal) THEN
              CALL smumps_sol_x_elt(mtype, id%N,
     &           id%NELT, id%ELTPTR(1),
     &           id%LELTVAR, id%ELTVAR(1),
     &           id%KEEP8(30), 
     &           id%A_ELT(1), sumr, keep(1),keep8(1) )
             ELSE
              CALL smumps_sol_scalx_elt(mtype, id%N,
     &           id%NELT, id%ELTPTR(1),
     &           id%LELTVAR, id%ELTVAR(1),
     &           id%KEEP8(30), 
     &           id%A_ELT(1),
     &           sumr, keep(1),keep8(1), id%COLSCA(1))
             ENDIF
            ENDIF
          ENDIF
      ELSE
          ALLOCATE( sumr_loc( id%N ), stat =allocok )
          IF (allocok .GT.0 ) THEN
             id%INFO(1)=-13
             id%INFO(2)=id%N
             RETURN
          ENDIF
          IF ( i_am_slave .and.
     &           id%KEEP8(29) .NE. 0 ) THEN
           IF (.NOT.lscal) THEN
              CALL smumps_sol_x(id%A_loc(1),
     &          id%KEEP8(29), id%N,
     &          id%IRN_loc(1), id%JCN_loc(1), 
     &          sumr_loc, id%KEEP(1),id%KEEP8(1), 
     &          eff_size_schur, id%SYM_PERM(1) )
           ELSE
              CALL smumps_scal_x(id%A_loc(1),
     &          id%KEEP8(29), id%N,
     &          id%IRN_loc(1), id%JCN_loc(1), 
     &          sumr_loc, id%KEEP(1),id%KEEP8(1),
     &          id%COLSCA(1),
     &          eff_size_schur, id%SYM_PERM(1) )
           ENDIF
          ELSE
           sumr_loc = zero
          ENDIF
          IF ( id%MYID .eq. master ) THEN
              CALL mpi_reduce( sumr_loc, sumr,
     &        id%N, mpi_real,
     &        mpi_sum,master,id%COMM, ierr)
          ELSE
              CALL mpi_reduce( sumr_loc, dummy,
     &        id%N, mpi_real,
     &        mpi_sum,master,id%COMM, ierr)
          END IF
        DEALLOCATE (sumr_loc)
      ENDIF
      IF ( id%MYID .eq. master ) THEN
        anorminf = real(zero)
        IF (lscal) THEN
         DO i = 1, id%N
          anorminf = max(abs(id%ROWSCA(i) * sumr(i)), 
     &                  anorminf)
         ENDDO
        ELSE
         DO i = 1, id%N
          anorminf = max(abs(sumr(i)), 
     &                  anorminf)
         ENDDO
        ENDIF
      ENDIF
      CALL mpi_bcast(anorminf, 1,
     &              mpi_real, master,
     &              id%COMM, ierr )
      IF (id%MYID .eq. master) DEALLOCATE (sumr)
      RETURN
      END SUBROUTINE smumps_anorminf