diff -Nur slicot/examples/TAB08ND.f slicot-new/examples/TAB08ND.f
--- slicot/examples/TAB08ND.f	2010-07-27 11:46:56.000000000 -0500
+++ slicot-new/examples/TAB08ND.f	2020-08-13 12:30:34.908489317 -0500
@@ -32,7 +32,7 @@
       INTEGER          INFZ(NMAX), IWORK(MPMAX+1), KRONL(NMAX+1),
      $                 KRONR(NMAX+1), LINFZ(NMAX)
 *     .. External Subroutines ..
-      EXTERNAL         AB08ND, DGEGV
+      EXTERNAL         AB08ND, DGGEV
 *     .. Intrinsic Functions ..
       INTRINSIC        MAX
 *     .. Executable Statements ..
@@ -114,7 +114,7 @@
 *                    Compute the invariant zeros of the given system.
 *                    Workspace: need 8*NU.
                      WRITE ( NOUT, FMT = 99983 )
-                     CALL DGEGV( 'No vectors', 'No vectors', NU, AF,
+                     CALL DGGEV( 'No vectors', 'No vectors', NU, AF,
      $                           LDAF, BF, LDBF, ALFR, ALFI, BETA, Q,
      $                           LDQ, Z, LDZ, DWORK, LDWORK, INFO )
 *
@@ -165,7 +165,7 @@
 *
 99999 FORMAT (' AB08ND EXAMPLE PROGRAM RESULTS',/1X)
 99998 FORMAT (' INFO on exit from AB08ND = ',I2)
-99997 FORMAT (' INFO on exit from DGEGV = ',I2)
+99997 FORMAT (' INFO on exit from DGGEV = ',I2)
 99994 FORMAT (' The left Kronecker indices of (A,C) are ',/(20(I3,2X)))
 99993 FORMAT (/' The system (A,C) is completely observable ')
 99992 FORMAT (/' The dimension of the observable subspace = ',I3)
diff -Nur slicot/examples/TAB08NZ.f slicot-new/examples/TAB08NZ.f
--- slicot/examples/TAB08NZ.f	2010-07-27 11:46:56.000000000 -0500
+++ slicot-new/examples/TAB08NZ.f	2020-08-13 23:36:42.686649320 -0500
@@ -38,7 +38,7 @@
       INTEGER          INFZ(NMAX), IWORK(MPMAX+1), KRONL(NMAX+1),
      $                 KRONR(NMAX+1), LINFZ(NMAX)
 *     .. External Subroutines ..
-      EXTERNAL         AB08NZ, ZGEGV
+      EXTERNAL         AB08NZ, ZGGEV
 *     .. Intrinsic Functions ..
       INTRINSIC        MAX
 *     .. Executable Statements ..
@@ -121,7 +121,7 @@
 *                    Complex Workspace: need 2*NU.
 *                    Real Workspace:    need 8*NU.
                      WRITE ( NOUT, FMT = 99983 )
-                     CALL ZGEGV( 'No vectors', 'No vectors', NU, AF,
+                     CALL ZGGEV( 'No vectors', 'No vectors', NU, AF,
      $                           LDAF, BF, LDBF, ALPHA, BETA, Q, LDQ,
      $                           Z, LDZ, ZWORK, LZWORK, DWORK, INFO )
 *
@@ -165,7 +165,7 @@
 *
 99999 FORMAT (' AB08NZ EXAMPLE PROGRAM RESULTS',/1X)
 99998 FORMAT (' INFO on exit from AB08NZ = ',I2)
-99997 FORMAT (' INFO on exit from ZGEGV = ',I2)
+99997 FORMAT (' INFO on exit from ZGGEV = ',I2)
 99994 FORMAT (' The left Kronecker indices of (A,C) are ',/(20(I3,2X)))
 99993 FORMAT (/' The system (A,C) is completely observable ')
 99992 FORMAT (/' The dimension of the observable subspace = ',I3)
diff -Nur slicot/examples/TAG08BD.f slicot-new/examples/TAG08BD.f
--- slicot/examples/TAG08BD.f	2010-07-27 11:46:56.000000000 -0500
+++ slicot-new/examples/TAG08BD.f	2020-08-13 21:45:45.526397715 -0500
@@ -31,7 +31,7 @@
      $                 IWORK(NMAX+MMAX), KRONL(LMAX+PMAX+1),
      $                 KRONR(NMAX+MMAX+1)
 *     .. External Subroutines ..
-      EXTERNAL         AG08BD, DGEGV, DLACPY
+      EXTERNAL         AG08BD, DGGEV, DLACPY
 *     .. Intrinsic Functions ..
       INTRINSIC        MAX
 *     .. Executable Statements ..
@@ -90,7 +90,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
    30                CONTINUE
-                     CALL DGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL DGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALFR, ALFI, BETA, Q,
      $                           LDQ, Z, LDZ, DWORK, LDWORK, INFO )
 *
@@ -146,7 +146,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
    70                CONTINUE
-                     CALL DGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL DGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALFR, ALFI, BETA, Q,
      $                           LDQ, Z, LDZ, DWORK, LDWORK, INFO )
 *
@@ -203,7 +203,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
   110                CONTINUE
-                     CALL DGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL DGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALFR, ALFI, BETA, Q,
      $                           LDQ, Z, LDZ, DWORK, LDWORK, INFO )
 *
@@ -258,7 +258,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
   140                CONTINUE
-                     CALL DGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL DGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALFR, ALFI, BETA, Q,
      $                           LDQ, Z, LDZ, DWORK, LDWORK, INFO )
 *
@@ -300,7 +300,7 @@
 *
 99999 FORMAT (' AG08BD EXAMPLE PROGRAM RESULTS',/1X)
 99998 FORMAT (' INFO on exit from AG08BD = ',I2)
-99997 FORMAT (' INFO on exit from DGEGV = ',I2)
+99997 FORMAT (' INFO on exit from DGGEV = ',I2)
 99996 FORMAT (/'Unobservable finite eigenvalues'/
      $         ' real  part     imag  part ')
 99995 FORMAT (/' The matrix Ef is ')
diff -Nur slicot/examples/TAG08BZ.f slicot-new/examples/TAG08BZ.f
--- slicot/examples/TAG08BZ.f	2010-07-27 11:46:56.000000000 -0500
+++ slicot-new/examples/TAG08BZ.f	2020-08-13 23:48:26.863342143 -0500
@@ -36,7 +36,7 @@
      $                 IWORK(NMAX+MMAX), KRONL(LMAX+PMAX+1),
      $                 KRONR(NMAX+MMAX+1)
 *     .. External Subroutines ..
-      EXTERNAL         AG08BZ, ZGEGV, ZLACPY
+      EXTERNAL         AG08BZ, ZGGEV, ZLACPY
 *     .. Intrinsic Functions ..
       INTRINSIC        MAX
 *     .. Executable Statements ..
@@ -96,7 +96,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
    30                CONTINUE
-                     CALL ZGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL ZGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALPHA, BETA, Q, LDQ,
      $                           Z, LDZ, ZWORK, LZWORK, DWORK, INFO )
 *
@@ -146,7 +146,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
    70                CONTINUE
-                     CALL ZGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL ZGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALPHA, BETA, Q, LDQ,
      $                           Z, LDZ, ZWORK, LZWORK, DWORK, INFO )
 *
@@ -197,7 +197,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
   110                CONTINUE
-                     CALL ZGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL ZGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALPHA, BETA, Q, LDQ,
      $                           Z, LDZ, ZWORK, LZWORK, DWORK, INFO )
 *
@@ -246,7 +246,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
   140                CONTINUE
-                     CALL ZGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL ZGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALPHA, BETA, Q, LDQ,
      $                           Z, LDZ, ZWORK, LZWORK, DWORK, INFO )
 *
@@ -281,7 +281,7 @@
 *
 99999 FORMAT (' AG08BZ EXAMPLE PROGRAM RESULTS',/1X)
 99998 FORMAT (' INFO on exit from AG08BZ = ',I2)
-99997 FORMAT (' INFO on exit from ZGEGV = ',I2)
+99997 FORMAT (' INFO on exit from ZGGEV = ',I2)
 99996 FORMAT (/' Unobservable finite eigenvalues'/
      $         ' real  part     imag  part ')
 99995 FORMAT (/' The matrix Ef is ')
diff -Nur slicot/examples77/TAB08ND.f slicot-new/examples77/TAB08ND.f
--- slicot/examples77/TAB08ND.f	2010-07-27 12:06:28.000000000 -0500
+++ slicot-new/examples77/TAB08ND.f	2020-08-13 23:59:53.073367645 -0500
@@ -32,7 +32,7 @@
       INTEGER          INFZ(NMAX), IWORK(MPMAX+1), KRONL(NMAX+1),
      $                 KRONR(NMAX+1), LINFZ(NMAX)
 *     .. External Subroutines ..
-      EXTERNAL         AB08ND, DGEGV
+      EXTERNAL         AB08ND, DGGEV
 *     .. Executable Statements ..
 *
       WRITE ( NOUT, FMT = 99999 )
@@ -112,7 +112,7 @@
 *                    Compute the invariant zeros of the given system.
 *                    Workspace: need 8*NU.
                      WRITE ( NOUT, FMT = 99983 )
-                     CALL DGEGV( 'No vectors', 'No vectors', NU, AF,
+                     CALL DGGEV( 'No vectors', 'No vectors', NU, AF,
      $                           LDAF, BF, LDBF, ALFR, ALFI, BETA, Q,
      $                           LDQ, Z, LDZ, DWORK, LDWORK, INFO )
 *
@@ -163,7 +163,7 @@
 *
 99999 FORMAT (' AB08ND EXAMPLE PROGRAM RESULTS',/1X)
 99998 FORMAT (' INFO on exit from AB08ND = ',I2)
-99997 FORMAT (' INFO on exit from DGEGV = ',I2)
+99997 FORMAT (' INFO on exit from DGGEV = ',I2)
 99994 FORMAT (' The left Kronecker indices of (A,C) are ',/(20(I3,2X)))
 99993 FORMAT (/' The system (A,C) is completely observable ')
 99992 FORMAT (/' The dimension of the observable subspace = ',I3)
diff -Nur slicot/examples77/TAB08NZ.f slicot-new/examples77/TAB08NZ.f
--- slicot/examples77/TAB08NZ.f	2010-07-27 12:06:28.000000000 -0500
+++ slicot-new/examples77/TAB08NZ.f	2020-08-14 00:04:37.360044912 -0500
@@ -32,7 +32,7 @@
       INTEGER          INFZ(NMAX), IWORK(MPMAX+1), KRONL(NMAX+1),
      $                 KRONR(NMAX+1), LINFZ(NMAX)
 *     .. External Subroutines ..
-      EXTERNAL         AB08NZ, ZGEGV
+      EXTERNAL         AB08NZ, ZGGEV
 *     .. Intrinsic Functions ..
       INTRINSIC        MAX
 *     .. Executable Statements ..
@@ -115,7 +115,7 @@
 *                    Complex Workspace: need 2*NU.
 *                    Real Workspace:    need 8*NU.
                      WRITE ( NOUT, FMT = 99983 )
-                     CALL ZGEGV( 'No vectors', 'No vectors', NU, AF,
+                     CALL ZGGEV( 'No vectors', 'No vectors', NU, AF,
      $                           LDAF, BF, LDBF, ALPHA, BETA, Q, LDQ,
      $                           Z, LDZ, ZWORK, LZWORK, DWORK, INFO )
 *
@@ -159,7 +159,7 @@
 *
 99999 FORMAT (' AB08NZ EXAMPLE PROGRAM RESULTS',/1X)
 99998 FORMAT (' INFO on exit from AB08NZ = ',I2)
-99997 FORMAT (' INFO on exit from ZGEGV = ',I2)
+99997 FORMAT (' INFO on exit from ZGGEV = ',I2)
 99994 FORMAT (' The left Kronecker indices of (A,C) are ',/(20(I3,2X)))
 99993 FORMAT (/' The system (A,C) is completely observable ')
 99992 FORMAT (/' The dimension of the observable subspace = ',I3)
diff -Nur slicot/examples77/TAG08BD.f slicot-new/examples77/TAG08BD.f
--- slicot/examples77/TAG08BD.f	2010-07-27 12:06:28.000000000 -0500
+++ slicot-new/examples77/TAG08BD.f	2020-08-14 00:02:27.393373386 -0500
@@ -31,7 +31,7 @@
      $                 IWORK(NMAX+MMAX), KRONL(LMAX+PMAX+1),
      $                 KRONR(NMAX+MMAX+1)
 *     .. External Subroutines ..
-      EXTERNAL         AG08BD, DGEGV, DLACPY
+      EXTERNAL         AG08BD, DGGEV, DLACPY
 *     .. Executable Statements ..
 *
       WRITE ( NOUT, FMT = 99999 )
@@ -88,7 +88,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
    30                CONTINUE
-                     CALL DGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL DGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALFR, ALFI, BETA, Q,
      $                           LDQ, Z, LDZ, DWORK, LDWORK, INFO )
 *
@@ -144,7 +144,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
    70                CONTINUE
-                     CALL DGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL DGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALFR, ALFI, BETA, Q,
      $                           LDQ, Z, LDZ, DWORK, LDWORK, INFO )
 *
@@ -201,7 +201,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
   110                CONTINUE
-                     CALL DGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL DGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALFR, ALFI, BETA, Q,
      $                           LDQ, Z, LDZ, DWORK, LDWORK, INFO )
 *
@@ -256,7 +256,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
   140                CONTINUE
-                     CALL DGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL DGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALFR, ALFI, BETA, Q,
      $                           LDQ, Z, LDZ, DWORK, LDWORK, INFO )
 *
@@ -298,7 +298,7 @@
 *
 99999 FORMAT (' AG08BD EXAMPLE PROGRAM RESULTS',/1X)
 99998 FORMAT (' INFO on exit from AG08BD = ',I2)
-99997 FORMAT (' INFO on exit from DGEGV = ',I2)
+99997 FORMAT (' INFO on exit from DGGEV = ',I2)
 99996 FORMAT (/'Unobservable finite eigenvalues'/
      $         ' real  part     imag  part ')
 99995 FORMAT (/' The matrix Ef is ')
diff -Nur slicot/examples77/TAG08BZ.f slicot-new/examples77/TAG08BZ.f
--- slicot/examples77/TAG08BZ.f	2010-07-27 12:06:28.000000000 -0500
+++ slicot-new/examples77/TAG08BZ.f	2020-08-14 00:06:30.673382447 -0500
@@ -34,7 +34,7 @@
      $                 IWORK(NMAX+MMAX), KRONL(LMAX+PMAX+1),
      $                 KRONR(NMAX+MMAX+1)
 *     .. External Subroutines ..
-      EXTERNAL         AG08BZ, ZGEGV, ZLACPY
+      EXTERNAL         AG08BZ, ZGGEV, ZLACPY
 *     .. Intrinsic Functions ..
       INTRINSIC        MAX
 *     .. Executable Statements ..
@@ -94,7 +94,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
    30                CONTINUE
-                     CALL ZGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL ZGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALPHA, BETA, Q, LDQ,
      $                           Z, LDZ, ZWORK, LZWORK, DWORK, INFO )
 *
@@ -144,7 +144,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
    70                CONTINUE
-                     CALL ZGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL ZGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALPHA, BETA, Q, LDQ,
      $                           Z, LDZ, ZWORK, LZWORK, DWORK, INFO )
 *
@@ -195,7 +195,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
   110                CONTINUE
-                     CALL ZGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL ZGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALPHA, BETA, Q, LDQ,
      $                           Z, LDZ, ZWORK, LZWORK, DWORK, INFO )
 *
@@ -244,7 +244,7 @@
                         WRITE ( NOUT, FMT = 99989 )
      $                        ( E(I,J), J = 1,NFZ )
   140                CONTINUE
-                     CALL ZGEGV( 'No vectors', 'No vectors', NFZ, A,
+                     CALL ZGGEV( 'No vectors', 'No vectors', NFZ, A,
      $                           LDA, E, LDE, ALPHA, BETA, Q, LDQ,
      $                           Z, LDZ, ZWORK, LZWORK, DWORK, INFO )
 *
@@ -279,7 +279,7 @@
 *
 99999 FORMAT (' AG08BZ EXAMPLE PROGRAM RESULTS',/1X)
 99998 FORMAT (' INFO on exit from AG08BZ = ',I2)
-99997 FORMAT (' INFO on exit from ZGEGV = ',I2)
+99997 FORMAT (' INFO on exit from ZGGEV = ',I2)
 99996 FORMAT (/' Unobservable finite eigenvalues'/
      $         ' real  part     imag  part ')
 99995 FORMAT (/' The matrix Ef is ')
diff -Nur slicot/src/AB08NX.f slicot-new/src/AB08NX.f
--- slicot/src/AB08NX.f	2010-07-27 11:14:47.000000000 -0500
+++ slicot-new/src/AB08NX.f	2020-08-13 13:20:19.665267163 -0500
@@ -199,7 +199,7 @@
       INTEGER           ILAENV
       EXTERNAL          ILAENV
 C     .. External Subroutines ..
-      EXTERNAL          DLAPMT, DLARFG, DLASET, DLATZM, DORMQR, DORMRQ,
+      EXTERNAL          DLAPMT, DTZRZF, DLASET, DORMRZ, DORMQR, DORMRQ,
      $                  MB03OY, MB03PY, XERBLA
 C     .. Intrinsic Functions ..
       INTRINSIC         INT, MAX, MIN
@@ -307,10 +307,12 @@
 C           Workspace: need   M+N-1.
 C
             DO 40 I1 = 1, SIGMA
-               CALL DLARFG( RO+1, ABCD(IROW,I1), ABCD(IROW+1,I1), 1, T )
-               CALL DLATZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1, T,
-     $                      ABCD(IROW,I1+1), ABCD(IROW+1,I1+1), LDABCD,
-     $                      DWORK )
+               CALL DTZRZF( 1, RO+1, ABCD(IROW,I1),  1, T,
+     $                      DWORK, LDWORK, INFO )
+               CALL DORMRZ( 'L', 'N', RO+1, MNU-I1, 1, RO,
+     $                      ABCD(IROW,I1), 1, T,
+     $                      ABCD(IROW,I1+1), LDABCD,
+     $                      DWORK , LDWORK, INFO )
                IROW = IROW + 1
    40       CONTINUE
             CALL DLASET( 'Lower', RO+SIGMA-1, SIGMA, ZERO, ZERO,
diff -Nur slicot/src/AB8NXZ.f slicot-new/src/AB8NXZ.f
--- slicot/src/AB8NXZ.f	2010-07-27 11:14:47.000000000 -0500
+++ slicot-new/src/AB8NXZ.f	2020-08-13 23:44:41.196667100 -0500
@@ -204,8 +204,8 @@
       INTEGER           ILAENV
       EXTERNAL          ILAENV
 C     .. External Subroutines ..
-      EXTERNAL          MB3OYZ, MB3PYZ, XERBLA, ZLAPMT, ZLARFG, ZLASET,
-     $                  ZLATZM, ZUNMQR, ZUNMRQ
+      EXTERNAL          MB3OYZ, MB3PYZ, XERBLA, ZLAPMT, ZTZRZF, ZLASET,
+     $                  ZUNMRZ, ZUNMQR, ZUNMRQ
 C     .. Intrinsic Functions ..
       INTRINSIC         DCONJG, INT, MAX, MIN
 C     .. Executable Statements ..
@@ -313,11 +313,12 @@
 C           CWorkspace: need   M+N-1.
 C
             DO 40 I1 = 1, SIGMA
-               CALL ZLARFG( RO+1, ABCD(IROW,I1), ABCD(IROW+1,I1), 1,
-     $                      TC )
-               CALL ZLATZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1,
-     $                      DCONJG( TC ), ABCD(IROW,I1+1),
-     $                      ABCD(IROW+1,I1+1), LDABCD, ZWORK )
+               CALL ZTZRZF( 1, RO+1, ABCD(IROW,I1), 1, TC, ZWORK,
+     $                      LZWORK, INFO )
+               CALL ZUNMRZ( 'L', 'N', RO+1, MNU-I1, 1, RO, 
+     $                      ABCD(IROW,I1), 1, DCONJG( TC ),
+     $                      ABCD(IROW,I1+1), LDABCD, ZWORK, LZWORK,
+     $                      INFO )
                IROW = IROW + 1
    40       CONTINUE
             CALL ZLASET( 'Lower', RO+SIGMA-1, SIGMA, ZERO, ZERO,
diff -Nur slicot/src/AG08BY.f slicot-new/src/AG08BY.f
--- slicot/src/AG08BY.f	2010-07-27 11:14:47.000000000 -0500
+++ slicot-new/src/AG08BY.f	2020-08-13 21:53:56.439749498 -0500
@@ -266,8 +266,8 @@
       DOUBLE PRECISION   DLAMCH, DNRM2
       EXTERNAL           DLAMCH, DNRM2, IDAMAX, ILAENV
 C     .. External Subroutines ..
-      EXTERNAL           DCOPY, DLAIC1, DLAPMT, DLARFG, DLARTG, DLASET,
-     $                   DLATZM, DORMQR, DROT, DSWAP, MB03OY, XERBLA
+      EXTERNAL           DCOPY, DLAIC1, DLAPMT, DTZRZF, DLARTG, DLASET,
+     $                   DORMRZ, DORMQR, DROT, DSWAP, MB03OY, XERBLA
 C     .. Intrinsic Functions ..
       INTRINSIC          ABS, DBLE, INT, MAX, MIN, SQRT
 C     .. Executable Statements ..
@@ -398,11 +398,12 @@
          IROW = NR
          DO 20 ICOL = 1, SIGMA
             IROW = IROW + 1
-            CALL DLARFG( RO+1, ABCD(IROW,ICOL), ABCD(IROW+1,ICOL), 1,
-     $                   T )
-            CALL DLATZM( 'L', RO+1, MNR-ICOL, ABCD(IROW+1,ICOL), 1, T,
-     $                   ABCD(IROW,ICOL+1), ABCD(IROW+1,ICOL+1),
-     $                   LDABCD, DWORK )
+            CALL DTZRZF( 1, RO+1, ABCD(IROW,ICOL), 1, T,
+     $                   DWORK, LDWORK, INFO )
+            CALL DORMRZ( 'L', 'N', RO+1, MNR-ICOL, 1, RO, 
+     $                   ABCD(IROW,ICOL), 1, T,
+     $                   ABCD(IROW,ICOL+1), LDABCD,
+     $                   DWORK, LDWORK, INFO )
             CALL DCOPY( PR-ICOL, DUM, 0, ABCD(IROW+1,ICOL), 1 )
    20    CONTINUE
          WRKOPT = MAX( WRKOPT, MN - 1 )
diff -Nur slicot/src/AG8BYZ.f slicot-new/src/AG8BYZ.f
--- slicot/src/AG8BYZ.f	2010-07-27 11:14:47.000000000 -0500
+++ slicot-new/src/AG8BYZ.f	2020-08-13 23:51:44.286682791 -0500
@@ -275,8 +275,8 @@
       DOUBLE PRECISION   DLAMCH, DZNRM2
       EXTERNAL           DLAMCH, DZNRM2, IDAMAX, ILAENV
 C     .. External Subroutines ..
-      EXTERNAL           MB3OYZ, XERBLA, ZCOPY, ZLAIC1, ZLAPMT, ZLARFG,
-     $                   ZLARTG, ZLASET, ZLATZM, ZROT, ZSWAP, ZUNMQR
+      EXTERNAL           MB3OYZ, XERBLA, ZCOPY, ZLAIC1, ZLAPMT, ZTZRZF,
+     $                   ZLARTG, ZLASET, ZUNMRZ, ZROT, ZSWAP, ZUNMQR
 C     .. Intrinsic Functions ..
       INTRINSIC          ABS, DBLE, DCONJG, INT, MAX, MIN, SQRT
 C     .. Executable Statements ..
@@ -407,11 +407,12 @@
          IROW = NR
          DO 20 ICOL = 1, SIGMA
             IROW = IROW + 1
-            CALL ZLARFG( RO+1, ABCD(IROW,ICOL), ABCD(IROW+1,ICOL), 1,
-     $                   TC )
-            CALL ZLATZM( 'L', RO+1, MNR-ICOL, ABCD(IROW+1,ICOL), 1,
-     $                   DCONJG( TC ), ABCD(IROW,ICOL+1),
-     $                   ABCD(IROW+1,ICOL+1), LDABCD, ZWORK )
+            CALL ZTZRZF( 1, RO+1, ABCD(IROW,ICOL), 1, TC,
+     $                   ZWORK, LZWORK, INFO )
+            CALL ZUNMRZ( 'L', 'N', RO+1, MNR-ICOL, 1, RO,
+     $                   ABCD(IROW,ICOL), 1, DCONJG( TC ),
+     $                   ABCD(IROW,ICOL+1), LDABCD,
+     $                   ZWORK, LZWORK, INFO )
             CALL ZCOPY( PR-ICOL, DUM, 0, ABCD(IROW+1,ICOL), 1 )
    20    CONTINUE
          WRKOPT = MAX( WRKOPT, MN - 1 )
diff -Nur slicot/src/SB01BY.f slicot-new/src/SB01BY.f
--- slicot/src/SB01BY.f	2010-07-27 11:14:48.000000000 -0500
+++ slicot-new/src/SB01BY.f	2020-08-13 21:36:33.463043580 -0500
@@ -113,16 +113,17 @@
 C     .. Array Arguments ..
       DOUBLE PRECISION  A(N,*), B(N,*), DWORK(*), F(M,*)
 C     .. Local Scalars ..
-      INTEGER           IR, J
+      INTEGER           IR, J, LDWORK
       DOUBLE PRECISION  ABSR, B1, B2, B21, C, C0, C1, C11, C12, C21,
      $                  C22, C3, C4, CS, CU, CV, DC0, DC2, DC3, DIFFR,
-     $                  R, RN, S12, S21, SIG, SN, SU, SV, TAU1, TAU2,
-     $                  WI, WI1, WR, WR1, X, Y, Z
+     $                  R, RN, S12, S21, SIG, SN, SU, SV, WI, WI1, WR,
+     $                  WR1, X, Y, Z
+      DOUBLE PRECISION  TAU1(N), TAU2(N)
 C     .. External Functions ..
       DOUBLE PRECISION  DLAMC3, DLAMCH
       EXTERNAL          DLAMC3, DLAMCH
 C     .. External Subroutines ..
-      EXTERNAL          DLANV2, DLARFG, DLASET, DLASV2, DLATZM, DROT
+      EXTERNAL          DLANV2, DTZRZF, DLASET, DLASV2, DORMRZ, DROT
 C     .. Intrinsic Functions ..
       INTRINSIC         ABS, MIN
 C     .. Executable Statements ..
@@ -130,12 +131,13 @@
 C     For efficiency reasons, the parameters are not checked.
 C
       INFO = 0
+      LDWORK = M 
       IF( N.EQ.1 ) THEN
 C
 C        The case N = 1.
 C
          IF( M.GT.1 )
-     $      CALL DLARFG( M, B(1,1), B(1,2), N, TAU1 )
+     $      CALL DTZRZF( 1, M, B(1,1), N, TAU1 , DWORK, LDWORK, INFO )
          B1 = B(1,1)
          IF( ABS( B1 ).LE.TOL ) THEN
 C
@@ -148,8 +150,9 @@
          F(1,1) = ( S - A(1,1) )/B1
          IF( M.GT.1 ) THEN
             CALL DLASET( 'Full', M-1, 1, ZERO, ZERO, F(2,1), M )
-            CALL DLATZM( 'Left', M, N, B(1,2), N, TAU1, F(1,1), F(2,1),
-     $                   M, DWORK )
+            CALL DORMRZ( 'Left', 'N', M, N, 1, M-1, 
+     $                   B(1,1), N, TAU1, F(1,1), M,
+     $                   DWORK, LDWORK, INFO )
          END IF
          RETURN
       END IF
@@ -184,13 +187,14 @@
 C        Postmultiply B with elementary Householder reflectors H1
 C        and H2.
 C
-         CALL DLARFG( M, B(1,1), B(1,2), N, TAU1 )
-         CALL DLATZM( 'Right', N-1, M, B(1,2), N, TAU1, B(2,1), B(2,2),
-     $                N, DWORK )
+         CALL DTZRZF( 1, M, B(1,1), N, TAU1, DWORK, LDWORK, INFO )
+         CALL DORMRZ( 'Right', 'N', N-1, M, N, M-1, 
+     $                B(1,1), N, TAU1, B(2,1), N,
+     $                DWORK, LDWORK, INFO )
          B1  = B(1,1)
          B21 = B(2,1)
          IF( M.GT.2 )
-     $      CALL DLARFG( M-1, B(2,2), B(2,3), N, TAU2 )
+     $      CALL DTZRZF( 1, M-1, B(2,2), N, TAU2, DWORK, LDWORK, INFO )
          B2  = B(2,2)
       END IF
 C
@@ -322,10 +326,12 @@
 C     Compute H1*H2*F.
 C
       IF( M.GT.2 )
-     $   CALL DLATZM( 'Left', M-1, N, B(2,3), N, TAU2, F(2,1), F(3,1),
-     $                M, DWORK )
-      CALL DLATZM( 'Left', M, N, B(1,2), N, TAU1, F(1,1), F(2,1), M,
-     $             DWORK )
+     $   CALL DORMRZ( 'Left', 'N', M-1, N, N, M-2,
+     $                B(2,2), N, TAU2, F(2,1), M,
+     $                DWORK, LDWORK, INFO )
+      CALL DORMRZ( 'Left', 'N', M, N, N, M-1,
+     $             B(1,1), N, TAU1, F(1,1), M,
+     $             DWORK, LDWORK, INFO )
 C
       RETURN
 C *** Last line of SB01BY ***
diff -Nur slicot/src/SB01FY.f slicot-new/src/SB01FY.f
--- slicot/src/SB01FY.f	2010-07-27 11:14:48.000000000 -0500
+++ slicot-new/src/SB01FY.f	2020-08-13 21:46:50.679733490 -0500
@@ -128,7 +128,7 @@
       DOUBLE PRECISION  DLAPY2, DLAPY3
       EXTERNAL          DLAPY2, DLAPY3
 C     .. External Subroutines ..
-      EXTERNAL          DLARFG, DLASET, DLATZM, DROTG, DTRTRI, MA02AD,
+      EXTERNAL          DTZRZF, DLASET, DORMRZ, DROTG, DTRTRI, MA02AD,
      $                  MB04OX, SB03OY
 C     .. Intrinsic Functions ..
       INTRINSIC         ABS, SQRT
@@ -179,9 +179,10 @@
 C        The N = 2 case.
 C
          IF( M.GT.1 ) THEN
-            CALL DLARFG( M, F(1,1), F(2,1), 1, TEMP )
-            CALL DLATZM( 'Left', M, N-1, F(2,1), 1, TEMP, F(1,2),
-     $                   F(2,2), LDF, V )
+            CALL DTZRZF( 1,  M, F(1,1), 1, TEMP, V, LDV, INFO )
+            CALL DORMRZ( 'Left', 'N', M, N-1, 1, M-1,
+     $                   F(1,1), 1, TEMP, F(1,2), LDF,
+     $                   V, LDV, INFO )
          END IF
          R11 = F(1,1)
          R12 = F(1,2)
diff -Nur slicot/src/SB04OD.f slicot-new/src/SB04OD.f
--- slicot/src/SB04OD.f	2010-07-27 11:14:48.000000000 -0500
+++ slicot-new/src/SB04OD.f	2020-08-13 22:50:18.553210749 -0500
@@ -413,7 +413,7 @@
       DOUBLE PRECISION  DLAMCH, DLANGE
       EXTERNAL          DLAMCH, DLANGE, LSAME
 C     .. External Subroutines ..
-      EXTERNAL          DCOPY, DGEGS, DGEMM, DGEMV, DLABAD, DLACPY,
+      EXTERNAL          DCOPY, DGGES, DGEMM, DGEMV, DLABAD, DLACPY,
      $                  DLASCL, DTGSYL, XERBLA
 C     .. Intrinsic Functions ..
       INTRINSIC         INT, MAX, SQRT
@@ -563,9 +563,10 @@
 C           Workspace:  need   7*M;
 C                       prefer 5*M + M*(NB+1).
 C
-            CALL DGEGS( 'Vectors left', 'Vectors right', M, A, LDA, D,
-     $                  LDD, DWORK, DWORK(M+1), DWORK(2*M+1), P, LDP, Q,
-     $                  LDQ, DWORK(3*M+1), LDWORK-3*M, INFO )
+            CALL DGGES( 'Vectors left', 'Vectors right', 'Not ordered', 
+     $                  NULL, M, A, LDA, D, LDD, NULL, DWORK,
+     $                  DWORK(M+1), DWORK(2*M+1), P, LDP, Q, LDQ,
+     $                  DWORK(3*M+1), LDWORK-3*M, NULL, INFO )
 C
 C           Undo scaling
 C
@@ -619,9 +620,10 @@
 C           Workspace:  need   7*N;
 C                       prefer 5*N + N*(NB+1).
 C
-            CALL DGEGS( 'Vectors left', 'Vectors right', N, B, LDB, E,
-     $                  LDE, DWORK, DWORK(N+1), DWORK(2*N+1), U, LDU, V,
-     $                  LDV, DWORK(3*N+1), LDWORK-3*N, INFO )
+            CALL DGGES( 'Vectors left', 'Vectors right', 'Not ordered',
+     $                  NULL, N, B, LDB, E, LDE, NULL, DWORK,
+     $                  DWORK(N+1), DWORK(2*N+1), U, LDU, V, LDV,
+     $                  DWORK(3*N+1), LDWORK-3*N, NULL, INFO )
 C
 C           Undo scaling
 C
diff -Nur slicot/src/SB06ND.f slicot-new/src/SB06ND.f
--- slicot/src/SB06ND.f	2010-07-27 11:14:48.000000000 -0500
+++ slicot-new/src/SB06ND.f	2020-08-20 16:15:51.248143212 -0500
@@ -181,13 +181,15 @@
       DOUBLE PRECISION  A(LDA,*), B(LDB,*), DWORK(*), F(LDF,*), U(LDU,*)
 C     .. Local Scalars ..
       INTEGER           J, J0, JCUR, JKCUR, JMKCUR, KCUR, KK, KMIN,
-     $                  KSTEP, MKCUR, NCONT
+     $                  KSTEP, MKCUR, NCONT, LDWORK
+      DOUBLE PRECISION  ATEMP(N), UTEMP(N), TAU
 C     .. External Subroutines ..
-      EXTERNAL          DCOPY, DGEMM, DLACPY, DLARFG, DLASET, DLATZM,
+      EXTERNAL          DCOPY, DGEMM, DLACPY, DLARFG, DLASET, DORMRZ,
      $                  DTRSM, XERBLA
 C     .. Executable Statements ..
 C
       INFO = 0
+      LDWORK = 2*N
 C
 C     Test the input scalar arguments.
 C
@@ -243,20 +245,29 @@
             DO 20 J = 1, KCUR
                JMKCUR = JCUR - KCUR
                CALL DCOPY( KCUR, A(JCUR,JMKCUR), LDA, F(1,JCUR), 1 )
-               CALL DLARFG( KCUR+1, A(JCUR,JCUR), F(1,JCUR), 1,
+               CALL DLARFG( KCUR+1, A(JCUR,JCUR), F(1,JCUR), 1, 
      $                      DWORK(JCUR) )
                CALL DLASET( 'Full', 1, KCUR, ZERO, ZERO, A(JCUR,JMKCUR),
      $                      LDA )
 C
 C              Backmultiply A and U with Ukk.
 C
-               CALL DLATZM( 'Right', JCUR-1, KCUR+1, F(1,JCUR), 1,
-     $                      DWORK(JCUR), A(1,JCUR), A(1,JMKCUR), LDA,
-     $                      DWORK )
-C
-               CALL DLATZM( 'Right', N, KCUR+1, F(1,JCUR), 1,
-     $                      DWORK(JCUR), U(1,JCUR), U(1,JMKCUR), LDU,
-     $                      DWORK(N+1) )
+               CALL DCOPY( JCUR-1, A(1,JCUR), 1, DWORK, 1 )
+               CALL DGEMV( 'N', JCUR-1, KCUR, ONE, A(1,JMKCUR), LDA, 
+     $                     F(1,JCUR), 1, ONE, DWORK, 1 )
+               CALL DAXPY( JCUR-1, -DWORK(JCUR), DWORK, 1,
+     $                     A(1,JCUR), 1 )
+               CALL DGER( JCUR-1, KCUR, -DWORK(JCUR), DWORK, 1,
+     $                    F(1,JCUR), 1, A(1,JMKCUR), LDA )
+C
+               CALL DCOPY( N, U(1,JCUR), 1, DWORK(N+1), 1 )
+               CALL DGEMV( 'N', N, KCUR, ONE, U(1,JMKCUR), LDA, 
+     $                     F(1,JCUR), 1, ONE, DWORK(N+1), 1 )
+               CALL DAXPY( N, -DWORK(JCUR), DWORK(N+1), 1,
+     $                     U(1,JCUR), 1 )
+               CALL DGER( N, KCUR, -DWORK(JCUR), DWORK(N+1),
+     $                    1, F(1,JCUR), 1, U(1,JMKCUR), LDA )
+
                JCUR = JCUR - 1
    20       CONTINUE
 C
@@ -291,9 +302,14 @@
                JCUR = JKCUR - KCUR
 C
                DO 60 J = 1, KCUR
-                  CALL DLATZM( 'Left', KCUR+1, N-JCUR+1, F(1,JKCUR), 1,
-     $                         DWORK(JKCUR), A(JKCUR,JCUR),
-     $                         A(JCUR,JCUR), LDA, DWORK(N+1) )
+                  CALL DCOPY( N-JCUR+1, A(JKCUR,JCUR), LDA,
+     $                        DWORK(N+1), 1)
+                  CALL DGEMV( 'T', KCUR, N-JCUR+1, ONE, A(JCUR,JCUR),
+     $                        LDA, F(1,JKCUR), 1, ONE, DWORK(N+1), 1 )
+                  CALL DAXPY( N-JCUR+1, -DWORK(JKCUR), DWORK(N+1), 1,
+     $                        A(JKCUR,JCUR), LDA ) 
+                  CALL DGER( KCUR, N-JCUR+1, -DWORK(JKCUR), F(1,JKCUR),
+     $                       1, DWORK(N+1), 1, A(JCUR,JCUR), LDA )
                   JCUR = JCUR - 1
                   JKCUR = JKCUR - 1
    60          CONTINUE
@@ -306,9 +322,13 @@
             JKCUR = JCUR + KCUR
 C
             DO 100 J = M, M - KCUR + 1, -1
-               CALL DLATZM( 'Left', KCUR+1, M-J+1, F(1,JKCUR), 1,
-     $                      DWORK(JKCUR), B(JKCUR,J), B(JCUR,J), LDB,
-     $                      DWORK(N+1) )
+               CALL DCOPY( M-J+1, B(JKCUR,J), LDB, DWORK(N+1), 1)
+               CALL DGEMV( 'T', KCUR, M-J+1, ONE, B(JCUR,J), LDB,
+     $                     F(1,JKCUR), 1, ONE, DWORK(N+1), 1 )
+               CALL DAXPY( M-J+1, -DWORK(JKCUR), DWORK(N+1), 1, 
+     $                     B(JKCUR,J), LDB ) 
+               CALL DGER( KCUR, M-J+1, -DWORK(JKCUR), F(1,JKCUR),
+     $                    1, DWORK(N+1), 1, B(JCUR,J), LDB )
                JCUR = JCUR - 1
                JKCUR = JKCUR - 1
   100       CONTINUE
diff -Nur slicot/src/SG03AD.f slicot-new/src/SG03AD.f
--- slicot/src/SG03AD.f	2010-07-27 11:14:48.000000000 -0500
+++ slicot-new/src/SG03AD.f	2020-08-13 22:56:31.469891587 -0500
@@ -214,7 +214,7 @@
 C                   quasitriangular form;
 C             = 2:  FACT = 'N' and the pencil A - lambda * E cannot be
 C                   reduced to generalized Schur form: LAPACK routine
-C                   DGEGS has failed to converge;
+C                   DGGES has failed to converge;
 C             = 3:  DICO = 'D' and the pencil A - lambda * E has a
 C                   pair of reciprocal eigenvalues. That is, lambda_i =
 C                   1/lambda_j for some i and j, where lambda_i and
@@ -392,7 +392,7 @@
       LOGICAL           LSAME
       EXTERNAL          DLAMCH, DNRM2, LSAME
 C     .. External Subroutines ..
-      EXTERNAL          DCOPY, DGEGS, DLACON, MB01RD, MB01RW, SG03AX,
+      EXTERNAL          DCOPY, DGGES, DLACON, MB01RD, MB01RW, SG03AX,
      $                  SG03AY, XERBLA
 C     .. Intrinsic Functions ..
       INTRINSIC         DBLE, INT, MAX, MIN
@@ -492,9 +492,9 @@
 C
 C        ( Workspace: >= MAX(1,4*N) )
 C
-         CALL DGEGS( 'Vectors', 'Vectors', N, A, LDA, E, LDE, ALPHAR,
-     $               ALPHAI, BETA, Q, LDQ, Z, LDZ, DWORK, LDWORK,
-     $               INFO1 )
+         CALL DGGES( 'Vectors', 'Vectors', 'N', NULL, N, A, LDA, E,
+     $               LDE, NULL, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ,
+     $               DWORK, LDWORK, NULL, INFO1 )
          IF ( INFO1 .NE. 0 ) THEN
             INFO = 2
             RETURN
diff -Nur slicot/src/SG03BD.f slicot-new/src/SG03BD.f
--- slicot/src/SG03BD.f	2010-07-27 11:14:48.000000000 -0500
+++ slicot-new/src/SG03BD.f	2020-08-13 22:58:24.363229236 -0500
@@ -203,7 +203,7 @@
 C                   eigenvalues are not conjugate complex;
 C             = 4:  FACT = 'N' and the pencil A - lambda * E cannot be
 C                   reduced to generalized Schur form: LAPACK routine
-C                   DGEGS has failed to converge;
+C                   DGGES has failed to converge;
 C             = 5:  DICO = 'C' and the pencil A - lambda * E is not
 C                   c-stable;
 C             = 6:  DICO = 'D' and the pencil A - lambda * E is not
@@ -472,7 +472,7 @@
       LOGICAL           LSAME
       EXTERNAL          DLAMCH, DLAPY2, LSAME
 C     .. External Subroutines ..
-      EXTERNAL          DCOPY, DGEGS, DGEMM, DGEMV, DGEQRF, DGERQF,
+      EXTERNAL          DCOPY, DGGES, DGEMM, DGEMV, DGEQRF, DGERQF,
      $                  DLACPY, DLAG2, DLASET, DSCAL, DTRMM, SG03BU,
      $                  SG03BV, XERBLA
 C     .. Intrinsic Functions ..
@@ -559,9 +559,9 @@
 C
 C        ( Workspace: >= MAX(1,4*N) )
 C
-         CALL DGEGS( 'Vectors', 'Vectors', N, A, LDA, E, LDE, ALPHAR,
-     $               ALPHAI, BETA, Q, LDQ, Z, LDZ, DWORK, LDWORK,
-     $               INFO1 )
+         CALL DGGES( 'Vectors', 'Vectors', 'N', NULL, N, A, LDA, E,
+     $               LDE, NULL, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ,
+     $               DWORK, LDWORK, NULL, INFO1 )
          IF ( INFO1 .NE. 0 ) THEN
             INFO = 4
             RETURN
diff -Nur slicot/src/TB01MD.f slicot-new/src/TB01MD.f
--- slicot/src/TB01MD.f	2010-07-27 11:14:48.000000000 -0500
+++ slicot-new/src/TB01MD.f	2020-08-13 23:12:53.886595888 -0500
@@ -182,18 +182,19 @@
 C     .. Local Scalars ..
       LOGICAL           LJOBA, LJOBI,  LUPLO
       INTEGER           II, J, M1, N1, NJ, PAR1, PAR2, PAR3, PAR4, PAR5,
-     $                  PAR6
+     $                  PAR6, LDWORK
       DOUBLE PRECISION  DZ
 C     .. External Functions ..
       LOGICAL           LSAME
       EXTERNAL          LSAME
 C     .. External Subroutines ..
-      EXTERNAL          DLARFG, DLASET, DLATZM, XERBLA
+      EXTERNAL          DTZRZF, DLASET, DORMRZ, XERBLA
 C     .. Intrinsic Functions ..
       INTRINSIC         MAX, MIN
 C     .. Executable Statements ..
 C
       INFO = 0
+      LDWORK = MAX(N,M-1)
       LUPLO = LSAME( UPLO, 'U' )
       LJOBI = LSAME( JOBU, 'I' )
       LJOBA = LJOBI.OR.LSAME( JOBU, 'U' )
@@ -258,29 +259,31 @@
             PAR5 = NJ
          END IF
 C
-         CALL DLARFG( NJ+1, B(PAR2,PAR1), B(PAR3,PAR1), 1, DZ )
+         CALL DTZRZF( 1, NJ+1, B(PAR2,PAR1), 1, DZ,
+     $                DWORK, LDWORK, INFO )
 C
 C        Update A.
 C
-         CALL DLATZM( 'Left', NJ+1, N, B(PAR3,PAR1), 1, DZ, A(PAR2,1),
-     $                A(PAR3,1), LDA, DWORK )
-         CALL DLATZM( 'Right', N, NJ+1, B(PAR3,PAR1), 1, DZ, A(1,PAR2),
-     $                A(1,PAR3), LDA, DWORK )
+         CALL DORMRZ( 'Left', 'N', NJ+1, N, 1, NJ, B(PAR2,PAR1), 1, DZ,
+     $                A(PAR2,1), LDA, DWORK, LDWORK, INFO )
+         CALL DORMRZ( 'Right', 'N', N, NJ+1, 1, NJ, B(PAR2,PAR1), 1, DZ,
+     $                A(1,PAR2), LDA, DWORK, LDWORK, INFO )
 C
          IF ( LJOBA ) THEN
 C
 C           Update U.
 C
-            CALL DLATZM( 'Right', N, NJ+1, B(PAR3,PAR1), 1, DZ,
-     $                   U(1,PAR2), U(1,PAR3), LDU, DWORK )
+            CALL DORMRZ( 'Right', 'N', N, NJ+1, 1, NJ, B(PAR2,PAR1),
+     $                   1, DZ, U(1,PAR2), LDU, DWORK, LDWORK, INFO )
          END IF
 C
          IF ( J.NE.M ) THEN
 C
 C           Update B
 C
-            CALL DLATZM( 'Left', NJ+1, PAR4-PAR3+1, B(PAR3,PAR1), 1, DZ,
-     $                   B(PAR2,PAR3), B(PAR3,PAR3), LDB, DWORK )
+            CALL DORMRZ( 'Left', 'N', NJ+1, PAR4-PAR3+1, 1, NJ, 
+     $                    B(PAR2,PAR1), 1, DZ, B(PAR2,PAR3), LDB,
+     $                    DWORK, LDWORK, INFO )
          END IF
 C
          DO 10 II = PAR3, PAR5
@@ -310,21 +313,22 @@
             PAR6 = N + M - J
          END IF
 C
-         CALL DLARFG( NJ+1, A(PAR2,PAR1), A(PAR3,PAR1), 1, DZ )
+         CALL DTZRZF( 1, NJ+1, A(PAR2,PAR1), 1, DZ, WORK, LDWORK, INFO )
 C
 C        Update A.
 C
-         CALL DLATZM( 'Left', NJ+1, PAR6-PAR5+1, A(PAR3,PAR1), 1, DZ,
-     $                A(PAR2,PAR5), A(PAR3,PAR5), LDA, DWORK )
-         CALL DLATZM( 'Right', N, NJ+1, A(PAR3,PAR1), 1, DZ,
-     $                A(1,PAR2), A(1,PAR3), LDA, DWORK )
+         CALL DORMRZ( 'Left', 'N', NJ+1, PAR6-PAR5+1, 1, NJ,
+     $                A(PAR2,PAR1), 1, DZ, A(PAR2,PAR5), LDA,
+     $                DWORK, LDWORK, INFO )
+         CALL DORMRZ( 'Right', 'N', N, NJ+1, 1, NJ, A(PAR2,PAR1), 1, DZ,
+     $                A(1,PAR2), LDA, DWORK, LDWORK, INFO )
 C
          IF ( LJOBA ) THEN
 C
 C           Update U.
 C
-            CALL DLATZM( 'Right', N, NJ+1, A(PAR3,PAR1), 1, DZ,
-     $                   U(1,PAR2), U(1,PAR3), LDU, DWORK )
+            CALL DORMRZ( 'Right', 'N', N, NJ+1, 1, NJ, A(PAR2,PAR1), 1,
+     $                   DZ, U(1,PAR2), LDU, DWORK, LDWORK, INFO )
          END IF
 C
          DO 30 II = PAR3, PAR4
diff -Nur slicot/src/TB01ND.f slicot-new/src/TB01ND.f
--- slicot/src/TB01ND.f	2010-07-27 11:14:48.000000000 -0500
+++ slicot-new/src/TB01ND.f	2020-08-20 16:15:43.534809593 -0500
@@ -189,18 +189,19 @@
 C     .. Local Scalars ..
       LOGICAL           LJOBA, LJOBI, LUPLO
       INTEGER           II, J, N1, NJ, P1, PAR1, PAR2, PAR3, PAR4, PAR5,
-     $                  PAR6
+     $                  PAR6, LDWORK
       DOUBLE PRECISION  DZ
 C     .. External Functions ..
       LOGICAL           LSAME
       EXTERNAL          LSAME
 C     .. External Subroutines ..
-      EXTERNAL          DLARFG, DLASET, DLATZM, XERBLA
+      EXTERNAL          DTZRZF, DLASET, DORMRZ, XERBLA
 C     .. Intrinsic Functions ..
       INTRINSIC         MAX, MIN
 C     .. Executable Statements ..
 C
       INFO = 0
+      LDWORK = MAX(N,P-1)
       LUPLO = LSAME( UPLO, 'U' )
       LJOBI = LSAME( JOBU, 'I' )
       LJOBA = LJOBI.OR.LSAME( JOBU, 'U' )
@@ -269,25 +270,41 @@
 C
 C        Update A.
 C
-         CALL DLATZM( 'Left', NJ+1, N, C(PAR1,PAR3), LDC, DZ, A(PAR2,1),
-     $                A(PAR3,1), LDA, DWORK )
-         CALL DLATZM( 'Right', N, NJ+1, C(PAR1,PAR3), LDC, DZ,
-     $                A(1,PAR2), A(1,PAR3), LDA, DWORK )
+         CALL DCOPY( N, A(PAR2,1), LDA, DWORK, 1)
+         CALL DGEMV( 'T', NJ, N, ONE, A(PAR3,1), LDA, C(PAR1,PAR3), LDC,
+     $               ONE, DWORK, 1 )
+         CALL DAXPY( N, -DZ, DWORK, 1, A(PAR2,1), LDA )
+         CALL DGER( NJ, N, -DZ, C(PAR1,PAR3), LDC, DWORK, 1,  
+     $              A(PAR3,1), LDA )
+         CALL DCOPY( N, A(1,PAR2), 1, DWORK, 1)
+         CALL DGEMV( 'N', N, NJ, ONE, A(1,PAR3), LDA, C(PAR1,PAR3), LDC,
+     $               ONE, DWORK, 1 )
+         CALL DAXPY( N, -DZ, DWORK, 1, A(1,PAR2), 1 )
+         CALL DGER( N, NJ, -DZ, DWORK, 1, C(PAR1,PAR3), LDC,
+     $              A(1,PAR3), LDA )
 C
          IF ( LJOBA ) THEN
 C
 C           Update U.
 C
-            CALL DLATZM( 'Right', N, NJ+1, C(PAR1,PAR3), LDC, DZ,
-     $                   U(1,PAR2), U(1,PAR3), LDU, DWORK )
+            CALL DCOPY( N, U(1,PAR2), 1, DWORK, 1)
+            CALL DGEMV( 'N', N, NJ, ONE, U(1,PAR3), LDU, C(PAR1,PAR3),
+     $                  LDC, ONE, DWORK, 1 )
+            CALL DAXPY( N, -DZ, DWORK, 1, U(1,PAR2), 1 )
+            CALL DGER( N, NJ, -DZ, DWORK , 1, C(PAR1,PAR3), LDC,
+     $              U(1,PAR3), LDU )
          END IF
 C
          IF ( J.NE.P ) THEN
 C
 C           Update C.
 C
-            CALL DLATZM( 'Right', PAR4-PAR3+1, NJ+1, C(PAR1,PAR3), LDC,
-     $                   DZ, C(PAR3,PAR2), C(PAR3,PAR3), LDC, DWORK )
+            CALL DCOPY( PAR4-PAR3+1, C(PAR3,PAR2), 1, DWORK, 1)
+            CALL DGEMV( 'N', PAR4-PAR3+1 , NJ, ONE, C(PAR3,PAR3), LDC,
+     $                  C(PAR1,PAR3), LDC, ONE, DWORK, 1 )
+            CALL DAXPY( PAR4-PAR3+1, -DZ, DWORK, 1, C(PAR3,PAR2), 1 )
+            CALL DGER( PAR4-PAR3+1, NJ, -DZ, DWORK , 1,
+     $                 C(PAR1,PAR3), LDC,C(PAR3,PAR3), LDU )
          END IF
 C
          DO 10 II = PAR3, PAR5
@@ -323,17 +340,32 @@
 C
 C           Update A.
 C
-            CALL DLATZM( 'Left', NJ+1, N, A(PAR1,PAR3), LDA, DZ,
-     $                   A(PAR2,1), A(PAR3,1), LDA, DWORK )
-            CALL DLATZM( 'Right', PAR6-PAR5+1, NJ+1, A(PAR1,PAR3), LDA,
-     $                    DZ, A(PAR5,PAR2), A(PAR5,PAR3), LDA, DWORK )
+            CALL DCOPY( N, A(PAR2,1), LDA, DWORK, 1)
+            CALL DGEMV( 'T', NJ, N, ONE, A(PAR3,1), LDA, A(PAR1,PAR3),
+     $                  LDA, ONE, DWORK, 1 )
+            CALL DAXPY( N, -DZ, DWORK, 1, A(PAR2,1), LDA )
+            CALL DGER( NJ, N, -DZ, A(PAR1,PAR3), LDA, DWORK, 1,
+     $                 A(PAR3,1), LDA )
+            CALL DCOPY( PAR6-PAR5+1, A(PAR5,PAR2), 1, DWORK, 1)
+            CALL DGEMV( 'N', PAR6-PAR5+1, NJ, ONE, A(PAR5,PAR3), LDA,
+     $                  A(PAR1,PAR3), LDA, ONE, DWORK, 1 )
+            CALL DAXPY( PAR6-PAR5+1, -DZ, DWORK, 1, A(PAR5,PAR2), 1 )
+            CALL DGER( PAR6-PAR5+1, NJ, -DZ, DWORK , 1, A(PAR1,PAR3),
+     $                 LDA, A(PAR5,PAR3), LDA )
 C
             IF ( LJOBA ) THEN
 C
 C              Update U.
 C
-               CALL DLATZM( 'Right', N, NJ+1, A(PAR1,PAR3), LDA, DZ,
-     $                      U(1,PAR2), U(1,PAR3), LDU, DWORK )
+               CALL DORMRZ( 'Right', 'N', N, NJ+1, 1, NJ,
+     $                      A(PAR1,PAR3), 1, DZ, U(1,PAR2), LDU,
+     $                      DWORK, LDWORK, INFO )
+               CALL DCOPY( N, U(1,PAR2), 1, DWORK, 1)
+               CALL DGEMV( 'N', N, NJ, ONE, U(1,PAR3), LDA,
+     $                     A(PAR1,PAR3), LDA, ONE, DWORK, 1 )
+               CALL DAXPY( N, -DZ, DWORK, 1, U(1,PAR2), 1 )
+               CALL DGER( N, NJ, -DZ, DWORK , 1, A(PAR1,PAR3),
+     $                 LDA, U(1,PAR3), LDA )
             END IF
 C
             DO 30 II = PAR3, PAR4
