Add level-3 BLAS triangular Sylvester equation solver

Author: angsch

SRC/CMakeLists.txt

@@ -306,7 +306,7 @@ set(DLASRC
    dlaqr0.f dlaqr1.f dlaqr2.f dlaqr3.f dlaqr4.f dlaqr5.f
    dlaqtr.f dlar1v.f dlar2v.f iladlr.f iladlc.f
    dlarf.f  dlarfb.f dlarfb_gett.f dlarfg.f dlarfgp.f dlarft.f dlarfx.f dlarfy.f
-   dlargv.f dlarrv.f dlartv.f
+   dlargv.f dlarmm.f dlarrv.f dlartv.f
    dlarz.f  dlarzb.f dlarzt.f dlaswp.f dlasy2.f
    dlasyf.f dlasyf_rook.f dlasyf_rk.f dlasyf_aa.f
    dlatbs.f dlatdf.f dlatps.f dlatrd.f dlatrs.f dlatrz.f dlauu2.f
@@ -342,7 +342,7 @@ set(DLASRC
    dtgsja.f dtgsna.f dtgsy2.f dtgsyl.f dtpcon.f dtprfs.f dtptri.f
    dtptrs.f
    dtrcon.f dtrevc.f dtrevc3.f dtrexc.f dtrrfs.f dtrsen.f dtrsna.f dtrsyl.f
-   dtrti2.f dtrtri.f dtrtrs.f dtzrzf.f dstemr.f
+   dtrsyl3.f dtrti2.f dtrtri.f dtrtrs.f dtzrzf.f dstemr.f
    dsgesv.f dsposv.f dlag2s.f slag2d.f dlat2s.f
    dlansf.f dpftrf.f dpftri.f dpftrs.f dsfrk.f dtfsm.f dtftri.f dtfttp.f
    dtfttr.f dtpttf.f dtpttr.f dtrttf.f dtrttp.f

SRC/Makefile

@@ -340,7 +340,7 @@ DLASRC = \
    dlaqr0.o dlaqr1.o dlaqr2.o dlaqr3.o dlaqr4.o dlaqr5.o \
    dlaqtr.o dlar1v.o dlar2v.o iladlr.o iladlc.o \
    dlarf.o  dlarfb.o dlarfb_gett.o dlarfg.o dlarfgp.o dlarft.o dlarfx.o dlarfy.o \
-   dlargv.o dlarrv.o dlartv.o \
+   dlargv.o dlarmm.o dlarrv.o dlartv.o \
    dlarz.o  dlarzb.o dlarzt.o dlaswp.o dlasy2.o \
    dlasyf.o dlasyf_rook.o dlasyf_rk.o \
    dlatbs.o dlatdf.o dlatps.o dlatrd.o dlatrs.o dlatrz.o dlauu2.o \
@@ -376,7 +376,7 @@ DLASRC = \
    dtgsja.o dtgsna.o dtgsy2.o dtgsyl.o dtpcon.o dtprfs.o dtptri.o \
    dtptrs.o \
    dtrcon.o dtrevc.o dtrevc3.o dtrexc.o dtrrfs.o dtrsen.o dtrsna.o dtrsyl.o \
-   dtrti2.o dtrtri.o dtrtrs.o dtzrzf.o dstemr.o \
+   dtrsyl3.o dtrti2.o dtrtri.o dtrtrs.o dtzrzf.o dstemr.o \
    dsgesv.o dsposv.o dlag2s.o slag2d.o dlat2s.o \
    dlansf.o dpftrf.o dpftri.o dpftrs.o dsfrk.o dtfsm.o dtftri.o dtfttp.o \
    dtfttr.o dtpttf.o dtpttr.o dtrttf.o dtrttp.o \

SRC/dlarmm.f

@@ -0,0 +1,103 @@
+*> \brief \b DLARMM 
+*
+* Definition:
+* ===========
+*
+*      DOUBLE PRECISION FUNCTION DLARMM( ANORM, BNORM, CNORM )
+*
+*     .. Scalar Arguments ..
+*      DOUBLE PRECISION   ANORM, BNORM, CNORM
+*     ..
+*
+*>  \par Purpose:
+*  =======
+*>
+*> \verbatim
+*>
+*> DLARMM returns a factor s in (0, 1] such that the linear updates
+*>
+*>    (s * C) - A * (s * B)  and  (s * C) - (s * A) * B
+*>
+*> cannot overflow, where A, B, and C are matrices of conforming
+*> dimensions.
+*>
+*> This is an auxiliary routine so there is no argument checking.
+*> \endverbatim
+*
+*  Arguments:
+*  =========
+*
+*> \param[in] ANORM
+*> \verbatim
+*>          ANORM is DOUBLE PRECISION
+*>          The infinity norm of A. ANORM >= 0.
+*>          The number of rows of the matrix A.  M >= 0.
+*> \endverbatim
+*>
+*> \param[in] BNORM
+*> \verbatim
+*>          BNORM is DOUBLE PRECISION
+*>          The infinity norm of B. BNORM >= 0.
+*> \endverbatim
+*>
+*> \param[in] CNORM
+*> \verbatim
+*>          CNORM is DOUBLE PRECISION
+*>          The infinity norm of C. CNORM >= 0.
+*> \endverbatim
+*>
+*>
+*  =====================================================================
+*>  Contributor:
+*>    Angelika Schwarz, Umea University, Sweden
+*>
+*  =====================================================================
+*>  References:
+*>    C. C. Kjelgaard Mikkelsen and L. Karlsson, Blocked Algorithms for
+*>    Robust Solution of Triangular Linear Systems. In: International
+*>    Conference on Parallel Processing and Applied Mathematics, pages
+*>    68--78. Springer, 2017.
+*>
+*> \ingroup OTHERauxiliary
+*  =====================================================================
+
+      DOUBLE PRECISION FUNCTION DLARMM( ANORM, BNORM, CNORM )
+      IMPLICIT NONE
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   ANORM, BNORM, CNORM
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, HALF
+      PARAMETER          ( ONE = 1.0D0, HALF = 0.5D+0 )
+*     ..
+*     .. Local Scalars ..
+       DOUBLE PRECISION   bignum, smlnum
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   dlamch
+      EXTERNAL           dlamch
+*     ..
+*     .. Executable Statements ..
+*
+*
+*     Determine machine dependent parameters to control overflow.
+*
+      SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' )
+      BIGNUM = ONE / SMLNUM
+*
+*     Compute a scale factor.
+*
+      DLARMM = ONE
+      IF( BNORM .LE. ONE ) THEN
+         IF( ANORM * BNORM .GT. BIGNUM - CNORM ) THEN
+            DLARMM = HALF
+         END IF
+      ELSE
+         IF( ANORM .GT. (BIGNUM - CNORM) / BNORM ) THEN
+            DLARMM = HALF / BNORM
+         END IF
+      END IF
+      RETURN
+*
+*     ==== End of DLARMM ====
+*
+      END