Add support for Sparse QR factorization

17 files changed, 6100 insertions, 1 deletions

diff --git a/lapack/CMakeLists.txt b/lapack/CMakeLists.txt
index 062845a3f..174e77d8d 100644
--- a/lapack/CMakeLists.txt
+++ b/lapack/CMakeLists.txt
@@ -18,6 +18,13 @@ single.cpp double.cpp complex_single.cpp complex_double.cpp ../blas/xerbla.cpp
 
 if(EIGEN_Fortran_COMPILER_WORKS)
 
+set(EigenLapack_SRCS 
+slarft.f dlarft.f clarft.f  zlarft.f
+slarfb.f dlarfb.f clarfb.f zlarfb.f
+slarfg.f dlarfg.f clarfg.f zlarfg.f
+slarf.f dlarf.f clarf.f zlarf.f
+)
+
 get_filename_component(eigen_full_path_to_reference_to_reference_lapack "./reference/" ABSOLUTE)
 if(EXISTS ${eigen_full_path_to_reference_to_reference_lapack})
 set(EigenLapack_SRCS ${EigenLapack_SRCS}
@@ -26,7 +33,6 @@ set(EigenLapack_SRCS ${EigenLapack_SRCS}
 # reference/dgetrf.f reference/cgetrf.f reference/sgetrf.f reference/zgetrf.f
 # reference/cgetrs.f reference/dgetrs.f reference/sgetrs.f reference/zgetrs.f
 # reference/dsyev.f  reference/ssyev.f
-
 reference/dlamch.f  reference/ilaver.f  reference/lsame.f  reference/slamch.f  reference/second_NONE.f  reference/dsecnd_NONE.f
 reference/cbdsqr.f               reference/ctbrfs.f               reference/dorml2.f               reference/sla_porfsx_extended.f  reference/zggglm.f
 reference/cgbbrd.f               reference/ctbtrs.f               reference/dormlq.f               reference/sla_porpvgrw.f         reference/zgghrd.f
diff --git a/lapack/clarf.f b/lapack/clarf.f
new file mode 100644
index 000000000..ca0328fb5
--- /dev/null
+++ b/lapack/clarf.f
@@ -0,0 +1,232 @@
+*> \brief \b CLARF
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download CLARF + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarf.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarf.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarf.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          SIDE
+*       INTEGER            INCV, LDC, M, N
+*       COMPLEX            TAU
+*       ..
+*       .. Array Arguments ..
+*       COMPLEX            C( LDC, * ), V( * ), WORK( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CLARF applies a complex elementary reflector H to a complex M-by-N
+*> matrix C, from either the left or the right. H is represented in the
+*> form
+*>
+*>       H = I - tau * v * v**H
+*>
+*> where tau is a complex scalar and v is a complex vector.
+*>
+*> If tau = 0, then H is taken to be the unit matrix.
+*>
+*> To apply H**H (the conjugate transpose of H), supply conjg(tau) instead
+*> tau.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*>          SIDE is CHARACTER*1
+*>          = 'L': form  H * C
+*>          = 'R': form  C * H
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*>          M is INTEGER
+*>          The number of rows of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of columns of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is COMPLEX array, dimension
+*>                     (1 + (M-1)*abs(INCV)) if SIDE = 'L'
+*>                  or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
+*>          The vector v in the representation of H. V is not used if
+*>          TAU = 0.
+*> \endverbatim
+*>
+*> \param[in] INCV
+*> \verbatim
+*>          INCV is INTEGER
+*>          The increment between elements of v. INCV <> 0.
+*> \endverbatim
+*>
+*> \param[in] TAU
+*> \verbatim
+*>          TAU is COMPLEX
+*>          The value tau in the representation of H.
+*> \endverbatim
+*>
+*> \param[in,out] C
+*> \verbatim
+*>          C is COMPLEX array, dimension (LDC,N)
+*>          On entry, the M-by-N matrix C.
+*>          On exit, C is overwritten by the matrix H * C if SIDE = 'L',
+*>          or C * H if SIDE = 'R'.
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX array, dimension
+*>                         (N) if SIDE = 'L'
+*>                      or (M) if SIDE = 'R'
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complexOTHERauxiliary
+*
+*  =====================================================================
+      SUBROUTINE CLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          SIDE
+      INTEGER            INCV, LDC, M, N
+      COMPLEX            TAU
+*     ..
+*     .. Array Arguments ..
+      COMPLEX            C( LDC, * ), V( * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX            ONE, ZERO
+      PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
+     $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            APPLYLEFT
+      INTEGER            I, LASTV, LASTC
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CGEMV, CGERC
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILACLR, ILACLC
+      EXTERNAL           LSAME, ILACLR, ILACLC
+*     ..
+*     .. Executable Statements ..
+*
+      APPLYLEFT = LSAME( SIDE, 'L' )
+      LASTV = 0
+      LASTC = 0
+      IF( TAU.NE.ZERO ) THEN
+!     Set up variables for scanning V.  LASTV begins pointing to the end
+!     of V.
+         IF( APPLYLEFT ) THEN
+            LASTV = M
+         ELSE
+            LASTV = N
+         END IF
+         IF( INCV.GT.0 ) THEN
+            I = 1 + (LASTV-1) * INCV
+         ELSE
+            I = 1
+         END IF
+!     Look for the last non-zero row in V.
+         DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO )
+            LASTV = LASTV - 1
+            I = I - INCV
+         END DO
+         IF( APPLYLEFT ) THEN
+!     Scan for the last non-zero column in C(1:lastv,:).
+            LASTC = ILACLC(LASTV, N, C, LDC)
+         ELSE
+!     Scan for the last non-zero row in C(:,1:lastv).
+            LASTC = ILACLR(M, LASTV, C, LDC)
+         END IF
+      END IF
+!     Note that lastc.eq.0 renders the BLAS operations null; no special
+!     case is needed at this level.
+      IF( APPLYLEFT ) THEN
+*
+*        Form  H * C
+*
+         IF( LASTV.GT.0 ) THEN
+*
+*           w(1:lastc,1) := C(1:lastv,1:lastc)**H * v(1:lastv,1)
+*
+            CALL CGEMV( 'Conjugate transpose', LASTV, LASTC, ONE,
+     $           C, LDC, V, INCV, ZERO, WORK, 1 )
+*
+*           C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)**H
+*
+            CALL CGERC( LASTV, LASTC, -TAU, V, INCV, WORK, 1, C, LDC )
+         END IF
+      ELSE
+*
+*        Form  C * H
+*
+         IF( LASTV.GT.0 ) THEN
+*
+*           w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1)
+*
+            CALL CGEMV( 'No transpose', LASTC, LASTV, ONE, C, LDC,
+     $           V, INCV, ZERO, WORK, 1 )
+*
+*           C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)**H
+*
+            CALL CGERC( LASTC, LASTV, -TAU, WORK, 1, V, INCV, C, LDC )
+         END IF
+      END IF
+      RETURN
+*
+*     End of CLARF
+*
+      END
diff --git a/lapack/clarfb.f b/lapack/clarfb.f
new file mode 100644
index 000000000..40bbdf487
--- /dev/null
+++ b/lapack/clarfb.f
@@ -0,0 +1,771 @@
+*> \brief \b CLARFB
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download CLARFB + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarfb.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarfb.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarfb.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
+*                          T, LDT, C, LDC, WORK, LDWORK )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          DIRECT, SIDE, STOREV, TRANS
+*       INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
+*       ..
+*       .. Array Arguments ..
+*       COMPLEX            C( LDC, * ), T( LDT, * ), V( LDV, * ),
+*      $                   WORK( LDWORK, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CLARFB applies a complex block reflector H or its transpose H**H to a
+*> complex M-by-N matrix C, from either the left or the right.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*>          SIDE is CHARACTER*1
+*>          = 'L': apply H or H**H from the Left
+*>          = 'R': apply H or H**H from the Right
+*> \endverbatim
+*>
+*> \param[in] TRANS
+*> \verbatim
+*>          TRANS is CHARACTER*1
+*>          = 'N': apply H (No transpose)
+*>          = 'C': apply H**H (Conjugate transpose)
+*> \endverbatim
+*>
+*> \param[in] DIRECT
+*> \verbatim
+*>          DIRECT is CHARACTER*1
+*>          Indicates how H is formed from a product of elementary
+*>          reflectors
+*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
+*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
+*> \endverbatim
+*>
+*> \param[in] STOREV
+*> \verbatim
+*>          STOREV is CHARACTER*1
+*>          Indicates how the vectors which define the elementary
+*>          reflectors are stored:
+*>          = 'C': Columnwise
+*>          = 'R': Rowwise
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*>          M is INTEGER
+*>          The number of rows of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of columns of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] K
+*> \verbatim
+*>          K is INTEGER
+*>          The order of the matrix T (= the number of elementary
+*>          reflectors whose product defines the block reflector).
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is COMPLEX array, dimension
+*>                                (LDV,K) if STOREV = 'C'
+*>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
+*>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
+*>          The matrix V. See Further Details.
+*> \endverbatim
+*>
+*> \param[in] LDV
+*> \verbatim
+*>          LDV is INTEGER
+*>          The leading dimension of the array V.
+*>          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
+*>          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
+*>          if STOREV = 'R', LDV >= K.
+*> \endverbatim
+*>
+*> \param[in] T
+*> \verbatim
+*>          T is COMPLEX array, dimension (LDT,K)
+*>          The triangular K-by-K matrix T in the representation of the
+*>          block reflector.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*>          LDT is INTEGER
+*>          The leading dimension of the array T. LDT >= K.
+*> \endverbatim
+*>
+*> \param[in,out] C
+*> \verbatim
+*>          C is COMPLEX array, dimension (LDC,N)
+*>          On entry, the M-by-N matrix C.
+*>          On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX array, dimension (LDWORK,K)
+*> \endverbatim
+*>
+*> \param[in] LDWORK
+*> \verbatim
+*>          LDWORK is INTEGER
+*>          The leading dimension of the array WORK.
+*>          If SIDE = 'L', LDWORK >= max(1,N);
+*>          if SIDE = 'R', LDWORK >= max(1,M).
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complexOTHERauxiliary
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  The shape of the matrix V and the storage of the vectors which define
+*>  the H(i) is best illustrated by the following example with n = 5 and
+*>  k = 3. The elements equal to 1 are not stored; the corresponding
+*>  array elements are modified but restored on exit. The rest of the
+*>  array is not used.
+*>
+*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
+*>
+*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
+*>                   ( v1  1    )                     (     1 v2 v2 v2 )
+*>                   ( v1 v2  1 )                     (        1 v3 v3 )
+*>                   ( v1 v2 v3 )
+*>                   ( v1 v2 v3 )
+*>
+*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
+*>
+*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
+*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
+*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
+*>                   (     1 v3 )
+*>                   (        1 )
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE CLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
+     $                   T, LDT, C, LDC, WORK, LDWORK )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          DIRECT, SIDE, STOREV, TRANS
+      INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
+*     ..
+*     .. Array Arguments ..
+      COMPLEX            C( LDC, * ), T( LDT, * ), V( LDV, * ),
+     $                   WORK( LDWORK, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX            ONE
+      PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
+*     ..
+*     .. Local Scalars ..
+      CHARACTER          TRANST
+      INTEGER            I, J, LASTV, LASTC
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILACLR, ILACLC
+      EXTERNAL           LSAME, ILACLR, ILACLC
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CCOPY, CGEMM, CLACGV, CTRMM
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          CONJG
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick return if possible
+*
+      IF( M.LE.0 .OR. N.LE.0 )
+     $   RETURN
+*
+      IF( LSAME( TRANS, 'N' ) ) THEN
+         TRANST = 'C'
+      ELSE
+         TRANST = 'N'
+      END IF
+*
+      IF( LSAME( STOREV, 'C' ) ) THEN
+*
+         IF( LSAME( DIRECT, 'F' ) ) THEN
+*
+*           Let  V =  ( V1 )    (first K rows)
+*                     ( V2 )
+*           where  V1  is unit lower triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**H * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILACLR( M, K, V, LDV ) )
+               LASTC = ILACLC( LASTV, N, C, LDC )
+*
+*              W := C**H * V  =  (C1**H * V1 + C2**H * V2)  (stored in WORK)
+*
+*              W := C1**H
+*
+               DO 10 J = 1, K
+                  CALL CCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
+                  CALL CLACGV( LASTC, WORK( 1, J ), 1 )
+   10          CONTINUE
+*
+*              W := W * V1
+*
+               CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2**H *V2
+*
+                  CALL CGEMM( 'Conjugate transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K, ONE, C( K+1, 1 ), LDC,
+     $                 V( K+1, 1 ), LDV, ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**H  or  W * T
+*
+               CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V * W**H
+*
+               IF( M.GT.K ) THEN
+*
+*                 C2 := C2 - V2 * W**H
+*
+                  CALL CGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTV-K, LASTC, K, -ONE, V( K+1, 1 ), LDV,
+     $                 WORK, LDWORK, ONE, C( K+1, 1 ), LDC )
+               END IF
+*
+*              W := W * V1**H
+*
+               CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W**H
+*
+               DO 30 J = 1, K
+                  DO 20 I = 1, LASTC
+                     C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) )
+   20             CONTINUE
+   30          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILACLR( N, K, V, LDV ) )
+               LASTC = ILACLR( M, LASTV, C, LDC )
+*
+*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
+*
+*              W := C1
+*
+               DO 40 J = 1, K
+                  CALL CCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
+   40          CONTINUE
+*
+*              W := W * V1
+*
+               CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2 * V2
+*
+                  CALL CGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**H
+*
+               CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V**H
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - W * V2**H
+*
+                  CALL CGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTC, LASTV-K, K,
+     $                 -ONE, WORK, LDWORK, V( K+1, 1 ), LDV,
+     $                 ONE, C( 1, K+1 ), LDC )
+               END IF
+*
+*              W := W * V1**H
+*
+               CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 60 J = 1, K
+                  DO 50 I = 1, LASTC
+                     C( I, J ) = C( I, J ) - WORK( I, J )
+   50             CONTINUE
+   60          CONTINUE
+            END IF
+*
+         ELSE
+*
+*           Let  V =  ( V1 )
+*                     ( V2 )    (last K rows)
+*           where  V2  is unit upper triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**H * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILACLR( M, K, V, LDV ) )
+               LASTC = ILACLC( LASTV, N, C, LDC )
+*
+*              W := C**H * V  =  (C1**H * V1 + C2**H * V2)  (stored in WORK)
+*
+*              W := C2**H
+*
+               DO 70 J = 1, K
+                  CALL CCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
+     $                 WORK( 1, J ), 1 )
+                  CALL CLACGV( LASTC, WORK( 1, J ), 1 )
+   70          CONTINUE
+*
+*              W := W * V2
+*
+               CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1**H*V1
+*
+                  CALL CGEMM( 'Conjugate transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**H  or  W * T
+*
+               CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V * W**H
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - V1 * W**H
+*
+                  CALL CGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2**H
+*
+               CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W**H
+*
+               DO 90 J = 1, K
+                  DO 80 I = 1, LASTC
+                     C( LASTV-K+J, I ) = C( LASTV-K+J, I ) -
+     $                               CONJG( WORK( I, J ) )
+   80             CONTINUE
+   90          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILACLR( N, K, V, LDV ) )
+               LASTC = ILACLR( M, LASTV, C, LDC )
+*
+*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
+*
+*              W := C2
+*
+               DO 100 J = 1, K
+                  CALL CCOPY( LASTC, C( 1, LASTV-K+J ), 1,
+     $                 WORK( 1, J ), 1 )
+  100          CONTINUE
+*
+*              W := W * V2
+*
+               CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1 * V1
+*
+                  CALL CGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**H
+*
+               CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V**H
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - W * V1**H
+*
+                  CALL CGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2**H
+*
+               CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W
+*
+               DO 120 J = 1, K
+                  DO 110 I = 1, LASTC
+                     C( I, LASTV-K+J ) = C( I, LASTV-K+J )
+     $                    - WORK( I, J )
+  110             CONTINUE
+  120          CONTINUE
+            END IF
+         END IF
+*
+      ELSE IF( LSAME( STOREV, 'R' ) ) THEN
+*
+         IF( LSAME( DIRECT, 'F' ) ) THEN
+*
+*           Let  V =  ( V1  V2 )    (V1: first K columns)
+*           where  V1  is unit upper triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**H * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILACLC( K, M, V, LDV ) )
+               LASTC = ILACLC( LASTV, N, C, LDC )
+*
+*              W := C**H * V**H  =  (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
+*
+*              W := C1**H
+*
+               DO 130 J = 1, K
+                  CALL CCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
+                  CALL CLACGV( LASTC, WORK( 1, J ), 1 )
+  130          CONTINUE
+*
+*              W := W * V1**H
+*
+               CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
+     $                     'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2**H*V2**H
+*
+                  CALL CGEMM( 'Conjugate transpose',
+     $                 'Conjugate transpose', LASTC, K, LASTV-K,
+     $                 ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**H  or  W * T
+*
+               CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V**H * W**H
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - V2**H * W**H
+*
+                  CALL CGEMM( 'Conjugate transpose',
+     $                 'Conjugate transpose', LASTV-K, LASTC, K,
+     $                 -ONE, V( 1, K+1 ), LDV, WORK, LDWORK,
+     $                 ONE, C( K+1, 1 ), LDC )
+               END IF
+*
+*              W := W * V1
+*
+               CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W**H
+*
+               DO 150 J = 1, K
+                  DO 140 I = 1, LASTC
+                     C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) )
+  140             CONTINUE
+  150          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILACLC( K, N, V, LDV ) )
+               LASTC = ILACLR( M, LASTV, C, LDC )
+*
+*              W := C * V**H  =  (C1*V1**H + C2*V2**H)  (stored in WORK)
+*
+*              W := C1
+*
+               DO 160 J = 1, K
+                  CALL CCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
+  160          CONTINUE
+*
+*              W := W * V1**H
+*
+               CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose',
+     $                     'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2 * V2**H
+*
+                  CALL CGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTC, K, LASTV-K, ONE, C( 1, K+1 ), LDC,
+     $                 V( 1, K+1 ), LDV, ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**H
+*
+               CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - W * V2
+*
+                  CALL CGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, LASTV-K, K,
+     $                 -ONE, WORK, LDWORK, V( 1, K+1 ), LDV,
+     $                 ONE, C( 1, K+1 ), LDC )
+               END IF
+*
+*              W := W * V1
+*
+               CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 180 J = 1, K
+                  DO 170 I = 1, LASTC
+                     C( I, J ) = C( I, J ) - WORK( I, J )
+  170             CONTINUE
+  180          CONTINUE
+*
+            END IF
+*
+         ELSE
+*
+*           Let  V =  ( V1  V2 )    (V2: last K columns)
+*           where  V2  is unit lower triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**H * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILACLC( K, M, V, LDV ) )
+               LASTC = ILACLC( LASTV, N, C, LDC )
+*
+*              W := C**H * V**H  =  (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
+*
+*              W := C2**H
+*
+               DO 190 J = 1, K
+                  CALL CCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
+     $                 WORK( 1, J ), 1 )
+                  CALL CLACGV( LASTC, WORK( 1, J ), 1 )
+  190          CONTINUE
+*
+*              W := W * V2**H
+*
+               CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1**H * V1**H
+*
+                  CALL CGEMM( 'Conjugate transpose',
+     $                 'Conjugate transpose', LASTC, K, LASTV-K,
+     $                 ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**H  or  W * T
+*
+               CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V**H * W**H
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - V1**H * W**H
+*
+                  CALL CGEMM( 'Conjugate transpose',
+     $                 'Conjugate transpose', LASTV-K, LASTC, K,
+     $                 -ONE, V, LDV, WORK, LDWORK, ONE, C, LDC )
+               END IF
+*
+*              W := W * V2
+*
+               CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W**H
+*
+               DO 210 J = 1, K
+                  DO 200 I = 1, LASTC
+                     C( LASTV-K+J, I ) = C( LASTV-K+J, I ) -
+     $                               CONJG( WORK( I, J ) )
+  200             CONTINUE
+  210          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILACLC( K, N, V, LDV ) )
+               LASTC = ILACLR( M, LASTV, C, LDC )
+*
+*              W := C * V**H  =  (C1*V1**H + C2*V2**H)  (stored in WORK)
+*
+*              W := C2
+*
+               DO 220 J = 1, K
+                  CALL CCOPY( LASTC, C( 1, LASTV-K+J ), 1,
+     $                 WORK( 1, J ), 1 )
+  220          CONTINUE
+*
+*              W := W * V2**H
+*
+               CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1 * V1**H
+*
+                  CALL CGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, ONE,
+     $                 WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**H
+*
+               CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - W * V1
+*
+                  CALL CGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2
+*
+               CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 240 J = 1, K
+                  DO 230 I = 1, LASTC
+                     C( I, LASTV-K+J ) = C( I, LASTV-K+J )
+     $                    - WORK( I, J )
+  230             CONTINUE
+  240          CONTINUE
+*
+            END IF
+*
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of CLARFB
+*
+      END
diff --git a/lapack/clarfg.f b/lapack/clarfg.f
new file mode 100644
index 000000000..d64f396c3
--- /dev/null
+++ b/lapack/clarfg.f
@@ -0,0 +1,203 @@
+*> \brief \b CLARFG
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download CLARFG + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarfg.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarfg.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarfg.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU )
+* 
+*       .. Scalar Arguments ..
+*       INTEGER            INCX, N
+*       COMPLEX            ALPHA, TAU
+*       ..
+*       .. Array Arguments ..
+*       COMPLEX            X( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CLARFG generates a complex elementary reflector H of order n, such
+*> that
+*>
+*>       H**H * ( alpha ) = ( beta ),   H**H * H = I.
+*>              (   x   )   (   0  )
+*>
+*> where alpha and beta are scalars, with beta real, and x is an
+*> (n-1)-element complex vector. H is represented in the form
+*>
+*>       H = I - tau * ( 1 ) * ( 1 v**H ) ,
+*>                     ( v )
+*>
+*> where tau is a complex scalar and v is a complex (n-1)-element
+*> vector. Note that H is not hermitian.
+*>
+*> If the elements of x are all zero and alpha is real, then tau = 0
+*> and H is taken to be the unit matrix.
+*>
+*> Otherwise  1 <= real(tau) <= 2  and  abs(tau-1) <= 1 .
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the elementary reflector.
+*> \endverbatim
+*>
+*> \param[in,out] ALPHA
+*> \verbatim
+*>          ALPHA is COMPLEX
+*>          On entry, the value alpha.
+*>          On exit, it is overwritten with the value beta.
+*> \endverbatim
+*>
+*> \param[in,out] X
+*> \verbatim
+*>          X is COMPLEX array, dimension
+*>                         (1+(N-2)*abs(INCX))
+*>          On entry, the vector x.
+*>          On exit, it is overwritten with the vector v.
+*> \endverbatim
+*>
+*> \param[in] INCX
+*> \verbatim
+*>          INCX is INTEGER
+*>          The increment between elements of X. INCX > 0.
+*> \endverbatim
+*>
+*> \param[out] TAU
+*> \verbatim
+*>          TAU is COMPLEX
+*>          The value tau.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complexOTHERauxiliary
+*
+*  =====================================================================
+      SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      INTEGER            INCX, N
+      COMPLEX            ALPHA, TAU
+*     ..
+*     .. Array Arguments ..
+      COMPLEX            X( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE, ZERO
+      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            J, KNT
+      REAL               ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM
+*     ..
+*     .. External Functions ..
+      REAL               SCNRM2, SLAMCH, SLAPY3
+      COMPLEX            CLADIV
+      EXTERNAL           SCNRM2, SLAMCH, SLAPY3, CLADIV
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, AIMAG, CMPLX, REAL, SIGN
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CSCAL, CSSCAL
+*     ..
+*     .. Executable Statements ..
+*
+      IF( N.LE.0 ) THEN
+         TAU = ZERO
+         RETURN
+      END IF
+*
+      XNORM = SCNRM2( N-1, X, INCX )
+      ALPHR = REAL( ALPHA )
+      ALPHI = AIMAG( ALPHA )
+*
+      IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN
+*
+*        H  =  I
+*
+         TAU = ZERO
+      ELSE
+*
+*        general case
+*
+         BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
+         SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' )
+         RSAFMN = ONE / SAFMIN
+*
+         KNT = 0
+         IF( ABS( BETA ).LT.SAFMIN ) THEN
+*
+*           XNORM, BETA may be inaccurate; scale X and recompute them
+*
+   10       CONTINUE
+            KNT = KNT + 1
+            CALL CSSCAL( N-1, RSAFMN, X, INCX )
+            BETA = BETA*RSAFMN
+            ALPHI = ALPHI*RSAFMN
+            ALPHR = ALPHR*RSAFMN
+            IF( ABS( BETA ).LT.SAFMIN )
+     $         GO TO 10
+*
+*           New BETA is at most 1, at least SAFMIN
+*
+            XNORM = SCNRM2( N-1, X, INCX )
+            ALPHA = CMPLX( ALPHR, ALPHI )
+            BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
+         END IF
+         TAU = CMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
+         ALPHA = CLADIV( CMPLX( ONE ), ALPHA-BETA )
+         CALL CSCAL( N-1, ALPHA, X, INCX )
+*
+*        If ALPHA is subnormal, it may lose relative accuracy
+*
+         DO 20 J = 1, KNT
+            BETA = BETA*SAFMIN
+ 20      CONTINUE
+         ALPHA = BETA
+      END IF
+*
+      RETURN
+*
+*     End of CLARFG
+*
+      END
diff --git a/lapack/clarft.f b/lapack/clarft.f
new file mode 100644
index 000000000..981447f77
--- /dev/null
+++ b/lapack/clarft.f
@@ -0,0 +1,328 @@
+*> \brief \b CLARFT
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download CLARFT + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarft.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarft.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarft.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          DIRECT, STOREV
+*       INTEGER            K, LDT, LDV, N
+*       ..
+*       .. Array Arguments ..
+*       COMPLEX            T( LDT, * ), TAU( * ), V( LDV, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> CLARFT forms the triangular factor T of a complex block reflector H
+*> of order n, which is defined as a product of k elementary reflectors.
+*>
+*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
+*>
+*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
+*>
+*> If STOREV = 'C', the vector which defines the elementary reflector
+*> H(i) is stored in the i-th column of the array V, and
+*>
+*>    H  =  I - V * T * V**H
+*>
+*> If STOREV = 'R', the vector which defines the elementary reflector
+*> H(i) is stored in the i-th row of the array V, and
+*>
+*>    H  =  I - V**H * T * V
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] DIRECT
+*> \verbatim
+*>          DIRECT is CHARACTER*1
+*>          Specifies the order in which the elementary reflectors are
+*>          multiplied to form the block reflector:
+*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
+*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
+*> \endverbatim
+*>
+*> \param[in] STOREV
+*> \verbatim
+*>          STOREV is CHARACTER*1
+*>          Specifies how the vectors which define the elementary
+*>          reflectors are stored (see also Further Details):
+*>          = 'C': columnwise
+*>          = 'R': rowwise
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the block reflector H. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] K
+*> \verbatim
+*>          K is INTEGER
+*>          The order of the triangular factor T (= the number of
+*>          elementary reflectors). K >= 1.
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is COMPLEX array, dimension
+*>                               (LDV,K) if STOREV = 'C'
+*>                               (LDV,N) if STOREV = 'R'
+*>          The matrix V. See further details.
+*> \endverbatim
+*>
+*> \param[in] LDV
+*> \verbatim
+*>          LDV is INTEGER
+*>          The leading dimension of the array V.
+*>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
+*> \endverbatim
+*>
+*> \param[in] TAU
+*> \verbatim
+*>          TAU is COMPLEX array, dimension (K)
+*>          TAU(i) must contain the scalar factor of the elementary
+*>          reflector H(i).
+*> \endverbatim
+*>
+*> \param[out] T
+*> \verbatim
+*>          T is COMPLEX array, dimension (LDT,K)
+*>          The k by k triangular factor T of the block reflector.
+*>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
+*>          lower triangular. The rest of the array is not used.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*>          LDT is INTEGER
+*>          The leading dimension of the array T. LDT >= K.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date April 2012
+*
+*> \ingroup complexOTHERauxiliary
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  The shape of the matrix V and the storage of the vectors which define
+*>  the H(i) is best illustrated by the following example with n = 5 and
+*>  k = 3. The elements equal to 1 are not stored.
+*>
+*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
+*>
+*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
+*>                   ( v1  1    )                     (     1 v2 v2 v2 )
+*>                   ( v1 v2  1 )                     (        1 v3 v3 )
+*>                   ( v1 v2 v3 )
+*>                   ( v1 v2 v3 )
+*>
+*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
+*>
+*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
+*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
+*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
+*>                   (     1 v3 )
+*>                   (        1 )
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
+*
+*  -- LAPACK auxiliary routine (version 3.4.1) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     April 2012
+*
+*     .. Scalar Arguments ..
+      CHARACTER          DIRECT, STOREV
+      INTEGER            K, LDT, LDV, N
+*     ..
+*     .. Array Arguments ..
+      COMPLEX            T( LDT, * ), TAU( * ), V( LDV, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX            ONE, ZERO
+      PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ),
+     $                   ZERO = ( 0.0E+0, 0.0E+0 ) )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, J, PREVLASTV, LASTV
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           CGEMV, CLACGV, CTRMV
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      IF( LSAME( DIRECT, 'F' ) ) THEN
+         PREVLASTV = N
+         DO I = 1, K
+            PREVLASTV = MAX( PREVLASTV, I )
+            IF( TAU( I ).EQ.ZERO ) THEN
+*
+*              H(i)  =  I
+*
+               DO J = 1, I
+                  T( J, I ) = ZERO
+               END DO
+            ELSE
+*
+*              general case
+*
+               IF( LSAME( STOREV, 'C' ) ) THEN
+*                 Skip any trailing zeros.
+                  DO LASTV = N, I+1, -1
+                     IF( V( LASTV, I ).NE.ZERO ) EXIT
+                  END DO
+                  DO J = 1, I-1
+                     T( J, I ) = -TAU( I ) * CONJG( V( I , J ) )
+                  END DO                     
+                  J = MIN( LASTV, PREVLASTV )
+*
+*                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i)
+*
+                  CALL CGEMV( 'Conjugate transpose', J-I, I-1,
+     $                        -TAU( I ), V( I+1, 1 ), LDV, 
+     $                        V( I+1, I ), 1,
+     $                        ONE, T( 1, I ), 1 )
+               ELSE
+*                 Skip any trailing zeros.
+                  DO LASTV = N, I+1, -1
+                     IF( V( I, LASTV ).NE.ZERO ) EXIT
+                  END DO
+                  DO J = 1, I-1
+                     T( J, I ) = -TAU( I ) * V( J , I )
+                  END DO                     
+                  J = MIN( LASTV, PREVLASTV )
+*
+*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H
+*
+                  CALL CGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ),
+     $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV,
+     $                        ONE, T( 1, I ), LDT )                  
+               END IF
+*
+*              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
+*
+               CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
+     $                     LDT, T( 1, I ), 1 )
+               T( I, I ) = TAU( I )
+               IF( I.GT.1 ) THEN
+                  PREVLASTV = MAX( PREVLASTV, LASTV )
+               ELSE
+                  PREVLASTV = LASTV
+               END IF
+            END IF
+         END DO
+      ELSE
+         PREVLASTV = 1
+         DO I = K, 1, -1
+            IF( TAU( I ).EQ.ZERO ) THEN
+*
+*              H(i)  =  I
+*
+               DO J = I, K
+                  T( J, I ) = ZERO
+               END DO
+            ELSE
+*
+*              general case
+*
+               IF( I.LT.K ) THEN
+                  IF( LSAME( STOREV, 'C' ) ) THEN
+*                    Skip any leading zeros.
+                     DO LASTV = 1, I-1
+                        IF( V( LASTV, I ).NE.ZERO ) EXIT
+                     END DO
+                     DO J = I+1, K
+                        T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) )
+                     END DO                        
+                     J = MAX( LASTV, PREVLASTV )
+*
+*                    T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i)
+*
+                     CALL CGEMV( 'Conjugate transpose', N-K+I-J, K-I,
+     $                           -TAU( I ), V( J, I+1 ), LDV, V( J, I ),
+     $                           1, ONE, T( I+1, I ), 1 )
+                  ELSE
+*                    Skip any leading zeros.
+                     DO LASTV = 1, I-1
+                        IF( V( I, LASTV ).NE.ZERO ) EXIT
+                     END DO
+                     DO J = I+1, K
+                        T( J, I ) = -TAU( I ) * V( J, N-K+I )
+                     END DO                      
+                     J = MAX( LASTV, PREVLASTV )
+*
+*                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H
+*
+                     CALL CGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ),
+     $                           V( I+1, J ), LDV, V( I, J ), LDV,
+     $                           ONE, T( I+1, I ), LDT )                     
+                  END IF
+*
+*                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
+*
+                  CALL CTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
+     $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
+                  IF( I.GT.1 ) THEN
+                     PREVLASTV = MIN( PREVLASTV, LASTV )
+                  ELSE
+                     PREVLASTV = LASTV
+                  END IF
+               END IF
+               T( I, I ) = TAU( I )
+            END IF
+         END DO
+      END IF
+      RETURN
+*
+*     End of CLARFT
+*
+      END
diff --git a/lapack/dlarf.f b/lapack/dlarf.f
new file mode 100644
index 000000000..2a82ff439
--- /dev/null
+++ b/lapack/dlarf.f
@@ -0,0 +1,227 @@
+*> \brief \b DLARF
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download DLARF + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarf.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarf.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarf.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          SIDE
+*       INTEGER            INCV, LDC, M, N
+*       DOUBLE PRECISION   TAU
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION   C( LDC, * ), V( * ), WORK( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DLARF applies a real elementary reflector H to a real m by n matrix
+*> C, from either the left or the right. H is represented in the form
+*>
+*>       H = I - tau * v * v**T
+*>
+*> where tau is a real scalar and v is a real vector.
+*>
+*> If tau = 0, then H is taken to be the unit matrix.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*>          SIDE is CHARACTER*1
+*>          = 'L': form  H * C
+*>          = 'R': form  C * H
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*>          M is INTEGER
+*>          The number of rows of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of columns of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is DOUBLE PRECISION array, dimension
+*>                     (1 + (M-1)*abs(INCV)) if SIDE = 'L'
+*>                  or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
+*>          The vector v in the representation of H. V is not used if
+*>          TAU = 0.
+*> \endverbatim
+*>
+*> \param[in] INCV
+*> \verbatim
+*>          INCV is INTEGER
+*>          The increment between elements of v. INCV <> 0.
+*> \endverbatim
+*>
+*> \param[in] TAU
+*> \verbatim
+*>          TAU is DOUBLE PRECISION
+*>          The value tau in the representation of H.
+*> \endverbatim
+*>
+*> \param[in,out] C
+*> \verbatim
+*>          C is DOUBLE PRECISION array, dimension (LDC,N)
+*>          On entry, the m by n matrix C.
+*>          On exit, C is overwritten by the matrix H * C if SIDE = 'L',
+*>          or C * H if SIDE = 'R'.
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is DOUBLE PRECISION array, dimension
+*>                         (N) if SIDE = 'L'
+*>                      or (M) if SIDE = 'R'
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup doubleOTHERauxiliary
+*
+*  =====================================================================
+      SUBROUTINE DLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          SIDE
+      INTEGER            INCV, LDC, M, N
+      DOUBLE PRECISION   TAU
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   C( LDC, * ), V( * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, ZERO
+      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            APPLYLEFT
+      INTEGER            I, LASTV, LASTC
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DGEMV, DGER
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILADLR, ILADLC
+      EXTERNAL           LSAME, ILADLR, ILADLC
+*     ..
+*     .. Executable Statements ..
+*
+      APPLYLEFT = LSAME( SIDE, 'L' )
+      LASTV = 0
+      LASTC = 0
+      IF( TAU.NE.ZERO ) THEN
+!     Set up variables for scanning V.  LASTV begins pointing to the end
+!     of V.
+         IF( APPLYLEFT ) THEN
+            LASTV = M
+         ELSE
+            LASTV = N
+         END IF
+         IF( INCV.GT.0 ) THEN
+            I = 1 + (LASTV-1) * INCV
+         ELSE
+            I = 1
+         END IF
+!     Look for the last non-zero row in V.
+         DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO )
+            LASTV = LASTV - 1
+            I = I - INCV
+         END DO
+         IF( APPLYLEFT ) THEN
+!     Scan for the last non-zero column in C(1:lastv,:).
+            LASTC = ILADLC(LASTV, N, C, LDC)
+         ELSE
+!     Scan for the last non-zero row in C(:,1:lastv).
+            LASTC = ILADLR(M, LASTV, C, LDC)
+         END IF
+      END IF
+!     Note that lastc.eq.0 renders the BLAS operations null; no special
+!     case is needed at this level.
+      IF( APPLYLEFT ) THEN
+*
+*        Form  H * C
+*
+         IF( LASTV.GT.0 ) THEN
+*
+*           w(1:lastc,1) := C(1:lastv,1:lastc)**T * v(1:lastv,1)
+*
+            CALL DGEMV( 'Transpose', LASTV, LASTC, ONE, C, LDC, V, INCV,
+     $           ZERO, WORK, 1 )
+*
+*           C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)**T
+*
+            CALL DGER( LASTV, LASTC, -TAU, V, INCV, WORK, 1, C, LDC )
+         END IF
+      ELSE
+*
+*        Form  C * H
+*
+         IF( LASTV.GT.0 ) THEN
+*
+*           w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1)
+*
+            CALL DGEMV( 'No transpose', LASTC, LASTV, ONE, C, LDC,
+     $           V, INCV, ZERO, WORK, 1 )
+*
+*           C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)**T
+*
+            CALL DGER( LASTC, LASTV, -TAU, WORK, 1, V, INCV, C, LDC )
+         END IF
+      END IF
+      RETURN
+*
+*     End of DLARF
+*
+      END
diff --git a/lapack/dlarfb.f b/lapack/dlarfb.f
new file mode 100644
index 000000000..206d3b268
--- /dev/null
+++ b/lapack/dlarfb.f
@@ -0,0 +1,762 @@
+*> \brief \b DLARFB
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download DLARFB + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfb.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfb.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfb.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
+*                          T, LDT, C, LDC, WORK, LDWORK )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          DIRECT, SIDE, STOREV, TRANS
+*       INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION   C( LDC, * ), T( LDT, * ), V( LDV, * ),
+*      $                   WORK( LDWORK, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DLARFB applies a real block reflector H or its transpose H**T to a
+*> real m by n matrix C, from either the left or the right.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*>          SIDE is CHARACTER*1
+*>          = 'L': apply H or H**T from the Left
+*>          = 'R': apply H or H**T from the Right
+*> \endverbatim
+*>
+*> \param[in] TRANS
+*> \verbatim
+*>          TRANS is CHARACTER*1
+*>          = 'N': apply H (No transpose)
+*>          = 'T': apply H**T (Transpose)
+*> \endverbatim
+*>
+*> \param[in] DIRECT
+*> \verbatim
+*>          DIRECT is CHARACTER*1
+*>          Indicates how H is formed from a product of elementary
+*>          reflectors
+*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
+*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
+*> \endverbatim
+*>
+*> \param[in] STOREV
+*> \verbatim
+*>          STOREV is CHARACTER*1
+*>          Indicates how the vectors which define the elementary
+*>          reflectors are stored:
+*>          = 'C': Columnwise
+*>          = 'R': Rowwise
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*>          M is INTEGER
+*>          The number of rows of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of columns of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] K
+*> \verbatim
+*>          K is INTEGER
+*>          The order of the matrix T (= the number of elementary
+*>          reflectors whose product defines the block reflector).
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is DOUBLE PRECISION array, dimension
+*>                                (LDV,K) if STOREV = 'C'
+*>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
+*>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
+*>          The matrix V. See Further Details.
+*> \endverbatim
+*>
+*> \param[in] LDV
+*> \verbatim
+*>          LDV is INTEGER
+*>          The leading dimension of the array V.
+*>          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
+*>          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
+*>          if STOREV = 'R', LDV >= K.
+*> \endverbatim
+*>
+*> \param[in] T
+*> \verbatim
+*>          T is DOUBLE PRECISION array, dimension (LDT,K)
+*>          The triangular k by k matrix T in the representation of the
+*>          block reflector.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*>          LDT is INTEGER
+*>          The leading dimension of the array T. LDT >= K.
+*> \endverbatim
+*>
+*> \param[in,out] C
+*> \verbatim
+*>          C is DOUBLE PRECISION array, dimension (LDC,N)
+*>          On entry, the m by n matrix C.
+*>          On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is DOUBLE PRECISION array, dimension (LDWORK,K)
+*> \endverbatim
+*>
+*> \param[in] LDWORK
+*> \verbatim
+*>          LDWORK is INTEGER
+*>          The leading dimension of the array WORK.
+*>          If SIDE = 'L', LDWORK >= max(1,N);
+*>          if SIDE = 'R', LDWORK >= max(1,M).
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup doubleOTHERauxiliary
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  The shape of the matrix V and the storage of the vectors which define
+*>  the H(i) is best illustrated by the following example with n = 5 and
+*>  k = 3. The elements equal to 1 are not stored; the corresponding
+*>  array elements are modified but restored on exit. The rest of the
+*>  array is not used.
+*>
+*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
+*>
+*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
+*>                   ( v1  1    )                     (     1 v2 v2 v2 )
+*>                   ( v1 v2  1 )                     (        1 v3 v3 )
+*>                   ( v1 v2 v3 )
+*>                   ( v1 v2 v3 )
+*>
+*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
+*>
+*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
+*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
+*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
+*>                   (     1 v3 )
+*>                   (        1 )
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
+     $                   T, LDT, C, LDC, WORK, LDWORK )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          DIRECT, SIDE, STOREV, TRANS
+      INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   C( LDC, * ), T( LDT, * ), V( LDV, * ),
+     $                   WORK( LDWORK, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE
+      PARAMETER          ( ONE = 1.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      CHARACTER          TRANST
+      INTEGER            I, J, LASTV, LASTC
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILADLR, ILADLC
+      EXTERNAL           LSAME, ILADLR, ILADLC
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DCOPY, DGEMM, DTRMM
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick return if possible
+*
+      IF( M.LE.0 .OR. N.LE.0 )
+     $   RETURN
+*
+      IF( LSAME( TRANS, 'N' ) ) THEN
+         TRANST = 'T'
+      ELSE
+         TRANST = 'N'
+      END IF
+*
+      IF( LSAME( STOREV, 'C' ) ) THEN
+*
+         IF( LSAME( DIRECT, 'F' ) ) THEN
+*
+*           Let  V =  ( V1 )    (first K rows)
+*                     ( V2 )
+*           where  V1  is unit lower triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**T * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILADLR( M, K, V, LDV ) )
+               LASTC = ILADLC( LASTV, N, C, LDC )
+*
+*              W := C**T * V  =  (C1**T * V1 + C2**T * V2)  (stored in WORK)
+*
+*              W := C1**T
+*
+               DO 10 J = 1, K
+                  CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
+   10          CONTINUE
+*
+*              W := W * V1
+*
+               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2**T *V2
+*
+                  CALL DGEMM( 'Transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**T  or  W * T
+*
+               CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V * W**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - V2 * W**T
+*
+                  CALL DGEMM( 'No transpose', 'Transpose',
+     $                 LASTV-K, LASTC, K,
+     $                 -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE,
+     $                 C( K+1, 1 ), LDC )
+               END IF
+*
+*              W := W * V1**T
+*
+               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W**T
+*
+               DO 30 J = 1, K
+                  DO 20 I = 1, LASTC
+                     C( J, I ) = C( J, I ) - WORK( I, J )
+   20             CONTINUE
+   30          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILADLR( N, K, V, LDV ) )
+               LASTC = ILADLR( M, LASTV, C, LDC )
+*
+*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
+*
+*              W := C1
+*
+               DO 40 J = 1, K
+                  CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
+   40          CONTINUE
+*
+*              W := W * V1
+*
+               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2 * V2
+*
+                  CALL DGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**T
+*
+               CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - W * V2**T
+*
+                  CALL DGEMM( 'No transpose', 'Transpose',
+     $                 LASTC, LASTV-K, K,
+     $                 -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE,
+     $                 C( 1, K+1 ), LDC )
+               END IF
+*
+*              W := W * V1**T
+*
+               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 60 J = 1, K
+                  DO 50 I = 1, LASTC
+                     C( I, J ) = C( I, J ) - WORK( I, J )
+   50             CONTINUE
+   60          CONTINUE
+            END IF
+*
+         ELSE
+*
+*           Let  V =  ( V1 )
+*                     ( V2 )    (last K rows)
+*           where  V2  is unit upper triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**T * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILADLR( M, K, V, LDV ) )
+               LASTC = ILADLC( LASTV, N, C, LDC )
+*
+*              W := C**T * V  =  (C1**T * V1 + C2**T * V2)  (stored in WORK)
+*
+*              W := C2**T
+*
+               DO 70 J = 1, K
+                  CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
+     $                 WORK( 1, J ), 1 )
+   70          CONTINUE
+*
+*              W := W * V2
+*
+               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1**T*V1
+*
+                  CALL DGEMM( 'Transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**T  or  W * T
+*
+               CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V * W**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - V1 * W**T
+*
+                  CALL DGEMM( 'No transpose', 'Transpose',
+     $                 LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2**T
+*
+               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W**T
+*
+               DO 90 J = 1, K
+                  DO 80 I = 1, LASTC
+                     C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J)
+   80             CONTINUE
+   90          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILADLR( N, K, V, LDV ) )
+               LASTC = ILADLR( M, LASTV, C, LDC )
+*
+*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
+*
+*              W := C2
+*
+               DO 100 J = 1, K
+                  CALL DCOPY( LASTC, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
+  100          CONTINUE
+*
+*              W := W * V2
+*
+               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1 * V1
+*
+                  CALL DGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**T
+*
+               CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - W * V1**T
+*
+                  CALL DGEMM( 'No transpose', 'Transpose',
+     $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2**T
+*
+               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W
+*
+               DO 120 J = 1, K
+                  DO 110 I = 1, LASTC
+                     C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J)
+  110             CONTINUE
+  120          CONTINUE
+            END IF
+         END IF
+*
+      ELSE IF( LSAME( STOREV, 'R' ) ) THEN
+*
+         IF( LSAME( DIRECT, 'F' ) ) THEN
+*
+*           Let  V =  ( V1  V2 )    (V1: first K columns)
+*           where  V1  is unit upper triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**T * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILADLC( K, M, V, LDV ) )
+               LASTC = ILADLC( LASTV, N, C, LDC )
+*
+*              W := C**T * V**T  =  (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
+*
+*              W := C1**T
+*
+               DO 130 J = 1, K
+                  CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
+  130          CONTINUE
+*
+*              W := W * V1**T
+*
+               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2**T*V2**T
+*
+                  CALL DGEMM( 'Transpose', 'Transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**T  or  W * T
+*
+               CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V**T * W**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - V2**T * W**T
+*
+                  CALL DGEMM( 'Transpose', 'Transpose',
+     $                 LASTV-K, LASTC, K,
+     $                 -ONE, V( 1, K+1 ), LDV, WORK, LDWORK,
+     $                 ONE, C( K+1, 1 ), LDC )
+               END IF
+*
+*              W := W * V1
+*
+               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W**T
+*
+               DO 150 J = 1, K
+                  DO 140 I = 1, LASTC
+                     C( J, I ) = C( J, I ) - WORK( I, J )
+  140             CONTINUE
+  150          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILADLC( K, N, V, LDV ) )
+               LASTC = ILADLR( M, LASTV, C, LDC )
+*
+*              W := C * V**T  =  (C1*V1**T + C2*V2**T)  (stored in WORK)
+*
+*              W := C1
+*
+               DO 160 J = 1, K
+                  CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
+  160          CONTINUE
+*
+*              W := W * V1**T
+*
+               CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2 * V2**T
+*
+                  CALL DGEMM( 'No transpose', 'Transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**T
+*
+               CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - W * V2
+*
+                  CALL DGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, LASTV-K, K,
+     $                 -ONE, WORK, LDWORK, V( 1, K+1 ), LDV,
+     $                 ONE, C( 1, K+1 ), LDC )
+               END IF
+*
+*              W := W * V1
+*
+               CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 180 J = 1, K
+                  DO 170 I = 1, LASTC
+                     C( I, J ) = C( I, J ) - WORK( I, J )
+  170             CONTINUE
+  180          CONTINUE
+*
+            END IF
+*
+         ELSE
+*
+*           Let  V =  ( V1  V2 )    (V2: last K columns)
+*           where  V2  is unit lower triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**T * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILADLC( K, M, V, LDV ) )
+               LASTC = ILADLC( LASTV, N, C, LDC )
+*
+*              W := C**T * V**T  =  (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
+*
+*              W := C2**T
+*
+               DO 190 J = 1, K
+                  CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
+     $                 WORK( 1, J ), 1 )
+  190          CONTINUE
+*
+*              W := W * V2**T
+*
+               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1**T * V1**T
+*
+                  CALL DGEMM( 'Transpose', 'Transpose',
+     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**T  or  W * T
+*
+               CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V**T * W**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - V1**T * W**T
+*
+                  CALL DGEMM( 'Transpose', 'Transpose',
+     $                 LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2
+*
+               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W**T
+*
+               DO 210 J = 1, K
+                  DO 200 I = 1, LASTC
+                     C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J)
+  200             CONTINUE
+  210          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILADLC( K, N, V, LDV ) )
+               LASTC = ILADLR( M, LASTV, C, LDC )
+*
+*              W := C * V**T  =  (C1*V1**T + C2*V2**T)  (stored in WORK)
+*
+*              W := C2
+*
+               DO 220 J = 1, K
+                  CALL DCOPY( LASTC, C( 1, LASTV-K+J ), 1,
+     $                 WORK( 1, J ), 1 )
+  220          CONTINUE
+*
+*              W := W * V2**T
+*
+               CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1 * V1**T
+*
+                  CALL DGEMM( 'No transpose', 'Transpose',
+     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**T
+*
+               CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - W * V1
+*
+                  CALL DGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2
+*
+               CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 240 J = 1, K
+                  DO 230 I = 1, LASTC
+                     C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J)
+  230             CONTINUE
+  240          CONTINUE
+*
+            END IF
+*
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of DLARFB
+*
+      END
diff --git a/lapack/dlarfg.f b/lapack/dlarfg.f
new file mode 100644
index 000000000..458ad2e05
--- /dev/null
+++ b/lapack/dlarfg.f
@@ -0,0 +1,196 @@
+*> \brief \b DLARFG
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download DLARFG + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfg.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfg.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfg.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU )
+* 
+*       .. Scalar Arguments ..
+*       INTEGER            INCX, N
+*       DOUBLE PRECISION   ALPHA, TAU
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION   X( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DLARFG generates a real elementary reflector H of order n, such
+*> that
+*>
+*>       H * ( alpha ) = ( beta ),   H**T * H = I.
+*>           (   x   )   (   0  )
+*>
+*> where alpha and beta are scalars, and x is an (n-1)-element real
+*> vector. H is represented in the form
+*>
+*>       H = I - tau * ( 1 ) * ( 1 v**T ) ,
+*>                     ( v )
+*>
+*> where tau is a real scalar and v is a real (n-1)-element
+*> vector.
+*>
+*> If the elements of x are all zero, then tau = 0 and H is taken to be
+*> the unit matrix.
+*>
+*> Otherwise  1 <= tau <= 2.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the elementary reflector.
+*> \endverbatim
+*>
+*> \param[in,out] ALPHA
+*> \verbatim
+*>          ALPHA is DOUBLE PRECISION
+*>          On entry, the value alpha.
+*>          On exit, it is overwritten with the value beta.
+*> \endverbatim
+*>
+*> \param[in,out] X
+*> \verbatim
+*>          X is DOUBLE PRECISION array, dimension
+*>                         (1+(N-2)*abs(INCX))
+*>          On entry, the vector x.
+*>          On exit, it is overwritten with the vector v.
+*> \endverbatim
+*>
+*> \param[in] INCX
+*> \verbatim
+*>          INCX is INTEGER
+*>          The increment between elements of X. INCX > 0.
+*> \endverbatim
+*>
+*> \param[out] TAU
+*> \verbatim
+*>          TAU is DOUBLE PRECISION
+*>          The value tau.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup doubleOTHERauxiliary
+*
+*  =====================================================================
+      SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      INTEGER            INCX, N
+      DOUBLE PRECISION   ALPHA, TAU
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   X( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, ZERO
+      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            J, KNT
+      DOUBLE PRECISION   BETA, RSAFMN, SAFMIN, XNORM
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLAMCH, DLAPY2, DNRM2
+      EXTERNAL           DLAMCH, DLAPY2, DNRM2
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, SIGN
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DSCAL
+*     ..
+*     .. Executable Statements ..
+*
+      IF( N.LE.1 ) THEN
+         TAU = ZERO
+         RETURN
+      END IF
+*
+      XNORM = DNRM2( N-1, X, INCX )
+*
+      IF( XNORM.EQ.ZERO ) THEN
+*
+*        H  =  I
+*
+         TAU = ZERO
+      ELSE
+*
+*        general case
+*
+         BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
+         SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' )
+         KNT = 0
+         IF( ABS( BETA ).LT.SAFMIN ) THEN
+*
+*           XNORM, BETA may be inaccurate; scale X and recompute them
+*
+            RSAFMN = ONE / SAFMIN
+   10       CONTINUE
+            KNT = KNT + 1
+            CALL DSCAL( N-1, RSAFMN, X, INCX )
+            BETA = BETA*RSAFMN
+            ALPHA = ALPHA*RSAFMN
+            IF( ABS( BETA ).LT.SAFMIN )
+     $         GO TO 10
+*
+*           New BETA is at most 1, at least SAFMIN
+*
+            XNORM = DNRM2( N-1, X, INCX )
+            BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
+         END IF
+         TAU = ( BETA-ALPHA ) / BETA
+         CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
+*
+*        If ALPHA is subnormal, it may lose relative accuracy
+*
+         DO 20 J = 1, KNT
+            BETA = BETA*SAFMIN
+ 20      CONTINUE
+         ALPHA = BETA
+      END IF
+*
+      RETURN
+*
+*     End of DLARFG
+*
+      END
diff --git a/lapack/dlarft.f b/lapack/dlarft.f
new file mode 100644
index 000000000..4b7550403
--- /dev/null
+++ b/lapack/dlarft.f
@@ -0,0 +1,326 @@
+*> \brief \b DLARFT
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download DLARFT + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarft.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarft.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarft.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          DIRECT, STOREV
+*       INTEGER            K, LDT, LDV, N
+*       ..
+*       .. Array Arguments ..
+*       DOUBLE PRECISION   T( LDT, * ), TAU( * ), V( LDV, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> DLARFT forms the triangular factor T of a real block reflector H
+*> of order n, which is defined as a product of k elementary reflectors.
+*>
+*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
+*>
+*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
+*>
+*> If STOREV = 'C', the vector which defines the elementary reflector
+*> H(i) is stored in the i-th column of the array V, and
+*>
+*>    H  =  I - V * T * V**T
+*>
+*> If STOREV = 'R', the vector which defines the elementary reflector
+*> H(i) is stored in the i-th row of the array V, and
+*>
+*>    H  =  I - V**T * T * V
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] DIRECT
+*> \verbatim
+*>          DIRECT is CHARACTER*1
+*>          Specifies the order in which the elementary reflectors are
+*>          multiplied to form the block reflector:
+*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
+*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
+*> \endverbatim
+*>
+*> \param[in] STOREV
+*> \verbatim
+*>          STOREV is CHARACTER*1
+*>          Specifies how the vectors which define the elementary
+*>          reflectors are stored (see also Further Details):
+*>          = 'C': columnwise
+*>          = 'R': rowwise
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the block reflector H. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] K
+*> \verbatim
+*>          K is INTEGER
+*>          The order of the triangular factor T (= the number of
+*>          elementary reflectors). K >= 1.
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is DOUBLE PRECISION array, dimension
+*>                               (LDV,K) if STOREV = 'C'
+*>                               (LDV,N) if STOREV = 'R'
+*>          The matrix V. See further details.
+*> \endverbatim
+*>
+*> \param[in] LDV
+*> \verbatim
+*>          LDV is INTEGER
+*>          The leading dimension of the array V.
+*>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
+*> \endverbatim
+*>
+*> \param[in] TAU
+*> \verbatim
+*>          TAU is DOUBLE PRECISION array, dimension (K)
+*>          TAU(i) must contain the scalar factor of the elementary
+*>          reflector H(i).
+*> \endverbatim
+*>
+*> \param[out] T
+*> \verbatim
+*>          T is DOUBLE PRECISION array, dimension (LDT,K)
+*>          The k by k triangular factor T of the block reflector.
+*>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
+*>          lower triangular. The rest of the array is not used.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*>          LDT is INTEGER
+*>          The leading dimension of the array T. LDT >= K.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date April 2012
+*
+*> \ingroup doubleOTHERauxiliary
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  The shape of the matrix V and the storage of the vectors which define
+*>  the H(i) is best illustrated by the following example with n = 5 and
+*>  k = 3. The elements equal to 1 are not stored.
+*>
+*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
+*>
+*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
+*>                   ( v1  1    )                     (     1 v2 v2 v2 )
+*>                   ( v1 v2  1 )                     (        1 v3 v3 )
+*>                   ( v1 v2 v3 )
+*>                   ( v1 v2 v3 )
+*>
+*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
+*>
+*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
+*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
+*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
+*>                   (     1 v3 )
+*>                   (        1 )
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
+*
+*  -- LAPACK auxiliary routine (version 3.4.1) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     April 2012
+*
+*     .. Scalar Arguments ..
+      CHARACTER          DIRECT, STOREV
+      INTEGER            K, LDT, LDV, N
+*     ..
+*     .. Array Arguments ..
+      DOUBLE PRECISION   T( LDT, * ), TAU( * ), V( LDV, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, ZERO
+      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, J, PREVLASTV, LASTV
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           DGEMV, DTRMV
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      IF( LSAME( DIRECT, 'F' ) ) THEN
+         PREVLASTV = N
+         DO I = 1, K
+            PREVLASTV = MAX( I, PREVLASTV )
+            IF( TAU( I ).EQ.ZERO ) THEN
+*
+*              H(i)  =  I
+*
+               DO J = 1, I
+                  T( J, I ) = ZERO
+               END DO
+            ELSE
+*
+*              general case
+*
+               IF( LSAME( STOREV, 'C' ) ) THEN
+*                 Skip any trailing zeros.
+                  DO LASTV = N, I+1, -1
+                     IF( V( LASTV, I ).NE.ZERO ) EXIT
+                  END DO
+                  DO J = 1, I-1
+                     T( J, I ) = -TAU( I ) * V( I , J )
+                  END DO   
+                  J = MIN( LASTV, PREVLASTV )
+*
+*                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i)
+*
+                  CALL DGEMV( 'Transpose', J-I, I-1, -TAU( I ), 
+     $                        V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE, 
+     $                        T( 1, I ), 1 )
+               ELSE
+*                 Skip any trailing zeros.
+                  DO LASTV = N, I+1, -1
+                     IF( V( I, LASTV ).NE.ZERO ) EXIT
+                  END DO
+                  DO J = 1, I-1
+                     T( J, I ) = -TAU( I ) * V( J , I )
+                  END DO   
+                  J = MIN( LASTV, PREVLASTV )
+*
+*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T
+*
+                  CALL DGEMV( 'No transpose', I-1, J-I, -TAU( I ),
+     $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV, ONE,
+     $                        T( 1, I ), 1 )
+               END IF
+*
+*              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
+*
+               CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
+     $                     LDT, T( 1, I ), 1 )
+               T( I, I ) = TAU( I )
+               IF( I.GT.1 ) THEN
+                  PREVLASTV = MAX( PREVLASTV, LASTV )
+               ELSE
+                  PREVLASTV = LASTV
+               END IF
+            END IF
+         END DO
+      ELSE
+         PREVLASTV = 1
+         DO I = K, 1, -1
+            IF( TAU( I ).EQ.ZERO ) THEN
+*
+*              H(i)  =  I
+*
+               DO J = I, K
+                  T( J, I ) = ZERO
+               END DO
+            ELSE
+*
+*              general case
+*
+               IF( I.LT.K ) THEN
+                  IF( LSAME( STOREV, 'C' ) ) THEN
+*                    Skip any leading zeros.
+                     DO LASTV = 1, I-1
+                        IF( V( LASTV, I ).NE.ZERO ) EXIT
+                     END DO
+                     DO J = I+1, K
+                        T( J, I ) = -TAU( I ) * V( N-K+I , J )
+                     END DO   
+                     J = MAX( LASTV, PREVLASTV )
+*
+*                    T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i)
+*
+                     CALL DGEMV( 'Transpose', N-K+I-J, K-I, -TAU( I ),
+     $                           V( J, I+1 ), LDV, V( J, I ), 1, ONE,
+     $                           T( I+1, I ), 1 )
+                  ELSE
+*                    Skip any leading zeros.
+                     DO LASTV = 1, I-1
+                        IF( V( I, LASTV ).NE.ZERO ) EXIT
+                     END DO
+                     DO J = I+1, K
+                        T( J, I ) = -TAU( I ) * V( J, N-K+I )
+                     END DO   
+                     J = MAX( LASTV, PREVLASTV )
+*
+*                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T
+*
+                     CALL DGEMV( 'No transpose', K-I, N-K+I-J,
+     $                    -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV,
+     $                    ONE, T( I+1, I ), 1 )
+                  END IF
+*
+*                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
+*
+                  CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
+     $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
+                  IF( I.GT.1 ) THEN
+                     PREVLASTV = MIN( PREVLASTV, LASTV )
+                  ELSE
+                     PREVLASTV = LASTV
+                  END IF
+               END IF
+               T( I, I ) = TAU( I )
+            END IF
+         END DO
+      END IF
+      RETURN
+*
+*     End of DLARFT
+*
+      END
diff --git a/lapack/slarf.f b/lapack/slarf.f
new file mode 100644
index 000000000..8a8ff308e
--- /dev/null
+++ b/lapack/slarf.f
@@ -0,0 +1,227 @@
+*> \brief \b SLARF
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download SLARF + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarf.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarf.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarf.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE SLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          SIDE
+*       INTEGER            INCV, LDC, M, N
+*       REAL               TAU
+*       ..
+*       .. Array Arguments ..
+*       REAL               C( LDC, * ), V( * ), WORK( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> SLARF applies a real elementary reflector H to a real m by n matrix
+*> C, from either the left or the right. H is represented in the form
+*>
+*>       H = I - tau * v * v**T
+*>
+*> where tau is a real scalar and v is a real vector.
+*>
+*> If tau = 0, then H is taken to be the unit matrix.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*>          SIDE is CHARACTER*1
+*>          = 'L': form  H * C
+*>          = 'R': form  C * H
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*>          M is INTEGER
+*>          The number of rows of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of columns of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is REAL array, dimension
+*>                     (1 + (M-1)*abs(INCV)) if SIDE = 'L'
+*>                  or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
+*>          The vector v in the representation of H. V is not used if
+*>          TAU = 0.
+*> \endverbatim
+*>
+*> \param[in] INCV
+*> \verbatim
+*>          INCV is INTEGER
+*>          The increment between elements of v. INCV <> 0.
+*> \endverbatim
+*>
+*> \param[in] TAU
+*> \verbatim
+*>          TAU is REAL
+*>          The value tau in the representation of H.
+*> \endverbatim
+*>
+*> \param[in,out] C
+*> \verbatim
+*>          C is REAL array, dimension (LDC,N)
+*>          On entry, the m by n matrix C.
+*>          On exit, C is overwritten by the matrix H * C if SIDE = 'L',
+*>          or C * H if SIDE = 'R'.
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is REAL array, dimension
+*>                         (N) if SIDE = 'L'
+*>                      or (M) if SIDE = 'R'
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup realOTHERauxiliary
+*
+*  =====================================================================
+      SUBROUTINE SLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          SIDE
+      INTEGER            INCV, LDC, M, N
+      REAL               TAU
+*     ..
+*     .. Array Arguments ..
+      REAL               C( LDC, * ), V( * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE, ZERO
+      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            APPLYLEFT
+      INTEGER            I, LASTV, LASTC
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SGEMV, SGER
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILASLR, ILASLC
+      EXTERNAL           LSAME, ILASLR, ILASLC
+*     ..
+*     .. Executable Statements ..
+*
+      APPLYLEFT = LSAME( SIDE, 'L' )
+      LASTV = 0
+      LASTC = 0
+      IF( TAU.NE.ZERO ) THEN
+!     Set up variables for scanning V.  LASTV begins pointing to the end
+!     of V.
+         IF( APPLYLEFT ) THEN
+            LASTV = M
+         ELSE
+            LASTV = N
+         END IF
+         IF( INCV.GT.0 ) THEN
+            I = 1 + (LASTV-1) * INCV
+         ELSE
+            I = 1
+         END IF
+!     Look for the last non-zero row in V.
+         DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO )
+            LASTV = LASTV - 1
+            I = I - INCV
+         END DO
+         IF( APPLYLEFT ) THEN
+!     Scan for the last non-zero column in C(1:lastv,:).
+            LASTC = ILASLC(LASTV, N, C, LDC)
+         ELSE
+!     Scan for the last non-zero row in C(:,1:lastv).
+            LASTC = ILASLR(M, LASTV, C, LDC)
+         END IF
+      END IF
+!     Note that lastc.eq.0 renders the BLAS operations null; no special
+!     case is needed at this level.
+      IF( APPLYLEFT ) THEN
+*
+*        Form  H * C
+*
+         IF( LASTV.GT.0 ) THEN
+*
+*           w(1:lastc,1) := C(1:lastv,1:lastc)**T * v(1:lastv,1)
+*
+            CALL SGEMV( 'Transpose', LASTV, LASTC, ONE, C, LDC, V, INCV,
+     $           ZERO, WORK, 1 )
+*
+*           C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)**T
+*
+            CALL SGER( LASTV, LASTC, -TAU, V, INCV, WORK, 1, C, LDC )
+         END IF
+      ELSE
+*
+*        Form  C * H
+*
+         IF( LASTV.GT.0 ) THEN
+*
+*           w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1)
+*
+            CALL SGEMV( 'No transpose', LASTC, LASTV, ONE, C, LDC,
+     $           V, INCV, ZERO, WORK, 1 )
+*
+*           C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)**T
+*
+            CALL SGER( LASTC, LASTV, -TAU, WORK, 1, V, INCV, C, LDC )
+         END IF
+      END IF
+      RETURN
+*
+*     End of SLARF
+*
+      END
diff --git a/lapack/slarfb.f b/lapack/slarfb.f
new file mode 100644
index 000000000..eb95990b3
--- /dev/null
+++ b/lapack/slarfb.f
@@ -0,0 +1,763 @@
+*> \brief \b SLARFB
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download SLARFB + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarfb.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarfb.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarfb.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE SLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
+*                          T, LDT, C, LDC, WORK, LDWORK )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          DIRECT, SIDE, STOREV, TRANS
+*       INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
+*       ..
+*       .. Array Arguments ..
+*       REAL               C( LDC, * ), T( LDT, * ), V( LDV, * ),
+*      $                   WORK( LDWORK, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> SLARFB applies a real block reflector H or its transpose H**T to a
+*> real m by n matrix C, from either the left or the right.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*>          SIDE is CHARACTER*1
+*>          = 'L': apply H or H**T from the Left
+*>          = 'R': apply H or H**T from the Right
+*> \endverbatim
+*>
+*> \param[in] TRANS
+*> \verbatim
+*>          TRANS is CHARACTER*1
+*>          = 'N': apply H (No transpose)
+*>          = 'T': apply H**T (Transpose)
+*> \endverbatim
+*>
+*> \param[in] DIRECT
+*> \verbatim
+*>          DIRECT is CHARACTER*1
+*>          Indicates how H is formed from a product of elementary
+*>          reflectors
+*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
+*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
+*> \endverbatim
+*>
+*> \param[in] STOREV
+*> \verbatim
+*>          STOREV is CHARACTER*1
+*>          Indicates how the vectors which define the elementary
+*>          reflectors are stored:
+*>          = 'C': Columnwise
+*>          = 'R': Rowwise
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*>          M is INTEGER
+*>          The number of rows of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of columns of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] K
+*> \verbatim
+*>          K is INTEGER
+*>          The order of the matrix T (= the number of elementary
+*>          reflectors whose product defines the block reflector).
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is REAL array, dimension
+*>                                (LDV,K) if STOREV = 'C'
+*>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
+*>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
+*>          The matrix V. See Further Details.
+*> \endverbatim
+*>
+*> \param[in] LDV
+*> \verbatim
+*>          LDV is INTEGER
+*>          The leading dimension of the array V.
+*>          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
+*>          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
+*>          if STOREV = 'R', LDV >= K.
+*> \endverbatim
+*>
+*> \param[in] T
+*> \verbatim
+*>          T is REAL array, dimension (LDT,K)
+*>          The triangular k by k matrix T in the representation of the
+*>          block reflector.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*>          LDT is INTEGER
+*>          The leading dimension of the array T. LDT >= K.
+*> \endverbatim
+*>
+*> \param[in,out] C
+*> \verbatim
+*>          C is REAL array, dimension (LDC,N)
+*>          On entry, the m by n matrix C.
+*>          On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T.
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is REAL array, dimension (LDWORK,K)
+*> \endverbatim
+*>
+*> \param[in] LDWORK
+*> \verbatim
+*>          LDWORK is INTEGER
+*>          The leading dimension of the array WORK.
+*>          If SIDE = 'L', LDWORK >= max(1,N);
+*>          if SIDE = 'R', LDWORK >= max(1,M).
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup realOTHERauxiliary
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  The shape of the matrix V and the storage of the vectors which define
+*>  the H(i) is best illustrated by the following example with n = 5 and
+*>  k = 3. The elements equal to 1 are not stored; the corresponding
+*>  array elements are modified but restored on exit. The rest of the
+*>  array is not used.
+*>
+*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
+*>
+*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
+*>                   ( v1  1    )                     (     1 v2 v2 v2 )
+*>                   ( v1 v2  1 )                     (        1 v3 v3 )
+*>                   ( v1 v2 v3 )
+*>                   ( v1 v2 v3 )
+*>
+*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
+*>
+*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
+*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
+*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
+*>                   (     1 v3 )
+*>                   (        1 )
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE SLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
+     $                   T, LDT, C, LDC, WORK, LDWORK )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          DIRECT, SIDE, STOREV, TRANS
+      INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
+*     ..
+*     .. Array Arguments ..
+      REAL               C( LDC, * ), T( LDT, * ), V( LDV, * ),
+     $                   WORK( LDWORK, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE
+      PARAMETER          ( ONE = 1.0E+0 )
+*     ..
+*     .. Local Scalars ..
+      CHARACTER          TRANST
+      INTEGER            I, J, LASTV, LASTC
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILASLR, ILASLC
+      EXTERNAL           LSAME, ILASLR, ILASLC
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SCOPY, SGEMM, STRMM
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick return if possible
+*
+      IF( M.LE.0 .OR. N.LE.0 )
+     $   RETURN
+*
+      IF( LSAME( TRANS, 'N' ) ) THEN
+         TRANST = 'T'
+      ELSE
+         TRANST = 'N'
+      END IF
+*
+      IF( LSAME( STOREV, 'C' ) ) THEN
+*
+         IF( LSAME( DIRECT, 'F' ) ) THEN
+*
+*           Let  V =  ( V1 )    (first K rows)
+*                     ( V2 )
+*           where  V1  is unit lower triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**T * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILASLR( M, K, V, LDV ) )
+               LASTC = ILASLC( LASTV, N, C, LDC )
+*
+*              W := C**T * V  =  (C1**T * V1 + C2**T * V2)  (stored in WORK)
+*
+*              W := C1**T
+*
+               DO 10 J = 1, K
+                  CALL SCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
+   10          CONTINUE
+*
+*              W := W * V1
+*
+               CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2**T *V2
+*
+                  CALL SGEMM( 'Transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**T  or  W * T
+*
+               CALL STRMM( 'Right', 'Upper', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V * W**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - V2 * W**T
+*
+                  CALL SGEMM( 'No transpose', 'Transpose',
+     $                 LASTV-K, LASTC, K,
+     $                 -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE,
+     $                 C( K+1, 1 ), LDC )
+               END IF
+*
+*              W := W * V1**T
+*
+               CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W**T
+*
+               DO 30 J = 1, K
+                  DO 20 I = 1, LASTC
+                     C( J, I ) = C( J, I ) - WORK( I, J )
+   20             CONTINUE
+   30          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILASLR( N, K, V, LDV ) )
+               LASTC = ILASLR( M, LASTV, C, LDC )
+*
+*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
+*
+*              W := C1
+*
+               DO 40 J = 1, K
+                  CALL SCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
+   40          CONTINUE
+*
+*              W := W * V1
+*
+               CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2 * V2
+*
+                  CALL SGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**T
+*
+               CALL STRMM( 'Right', 'Upper', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - W * V2**T
+*
+                  CALL SGEMM( 'No transpose', 'Transpose',
+     $                 LASTC, LASTV-K, K,
+     $                 -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE,
+     $                 C( 1, K+1 ), LDC )
+               END IF
+*
+*              W := W * V1**T
+*
+               CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 60 J = 1, K
+                  DO 50 I = 1, LASTC
+                     C( I, J ) = C( I, J ) - WORK( I, J )
+   50             CONTINUE
+   60          CONTINUE
+            END IF
+*
+         ELSE
+*
+*           Let  V =  ( V1 )
+*                     ( V2 )    (last K rows)
+*           where  V2  is unit upper triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**T * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILASLR( M, K, V, LDV ) )
+               LASTC = ILASLC( LASTV, N, C, LDC )
+*
+*              W := C**T * V  =  (C1**T * V1 + C2**T * V2)  (stored in WORK)
+*
+*              W := C2**T
+*
+               DO 70 J = 1, K
+                  CALL SCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
+     $                 WORK( 1, J ), 1 )
+   70          CONTINUE
+*
+*              W := W * V2
+*
+               CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1**T*V1
+*
+                  CALL SGEMM( 'Transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**T  or  W * T
+*
+               CALL STRMM( 'Right', 'Lower', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V * W**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - V1 * W**T
+*
+                  CALL SGEMM( 'No transpose', 'Transpose',
+     $                 LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2**T
+*
+               CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W**T
+*
+               DO 90 J = 1, K
+                  DO 80 I = 1, LASTC
+                     C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J)
+   80             CONTINUE
+   90          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILASLR( N, K, V, LDV ) )
+               LASTC = ILASLR( M, LASTV, C, LDC )
+*
+*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
+*
+*              W := C2
+*
+               DO 100 J = 1, K
+                  CALL SCOPY( LASTC, C( 1, N-K+J ), 1, WORK( 1, J ), 1 )
+  100          CONTINUE
+*
+*              W := W * V2
+*
+               CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1 * V1
+*
+                  CALL SGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**T
+*
+               CALL STRMM( 'Right', 'Lower', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - W * V1**T
+*
+                  CALL SGEMM( 'No transpose', 'Transpose',
+     $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2**T
+*
+               CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W
+*
+               DO 120 J = 1, K
+                  DO 110 I = 1, LASTC
+                     C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J)
+  110             CONTINUE
+  120          CONTINUE
+            END IF
+         END IF
+*
+      ELSE IF( LSAME( STOREV, 'R' ) ) THEN
+*
+         IF( LSAME( DIRECT, 'F' ) ) THEN
+*
+*           Let  V =  ( V1  V2 )    (V1: first K columns)
+*           where  V1  is unit upper triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**T * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILASLC( K, M, V, LDV ) )
+               LASTC = ILASLC( LASTV, N, C, LDC )
+*
+*              W := C**T * V**T  =  (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
+*
+*              W := C1**T
+*
+               DO 130 J = 1, K
+                  CALL SCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
+  130          CONTINUE
+*
+*              W := W * V1**T
+*
+               CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2**T*V2**T
+*
+                  CALL SGEMM( 'Transpose', 'Transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**T  or  W * T
+*
+               CALL STRMM( 'Right', 'Upper', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V**T * W**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - V2**T * W**T
+*
+                  CALL SGEMM( 'Transpose', 'Transpose',
+     $                 LASTV-K, LASTC, K,
+     $                 -ONE, V( 1, K+1 ), LDV, WORK, LDWORK,
+     $                 ONE, C( K+1, 1 ), LDC )
+               END IF
+*
+*              W := W * V1
+*
+               CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W**T
+*
+               DO 150 J = 1, K
+                  DO 140 I = 1, LASTC
+                     C( J, I ) = C( J, I ) - WORK( I, J )
+  140             CONTINUE
+  150          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILASLC( K, N, V, LDV ) )
+               LASTC = ILASLR( M, LASTV, C, LDC )
+*
+*              W := C * V**T  =  (C1*V1**T + C2*V2**T)  (stored in WORK)
+*
+*              W := C1
+*
+               DO 160 J = 1, K
+                  CALL SCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
+  160          CONTINUE
+*
+*              W := W * V1**T
+*
+               CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2 * V2**T
+*
+                  CALL SGEMM( 'No transpose', 'Transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**T
+*
+               CALL STRMM( 'Right', 'Upper', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - W * V2
+*
+                  CALL SGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, LASTV-K, K,
+     $                 -ONE, WORK, LDWORK, V( 1, K+1 ), LDV,
+     $                 ONE, C( 1, K+1 ), LDC )
+               END IF
+*
+*              W := W * V1
+*
+               CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 180 J = 1, K
+                  DO 170 I = 1, LASTC
+                     C( I, J ) = C( I, J ) - WORK( I, J )
+  170             CONTINUE
+  180          CONTINUE
+*
+            END IF
+*
+         ELSE
+*
+*           Let  V =  ( V1  V2 )    (V2: last K columns)
+*           where  V2  is unit lower triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**T * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILASLC( K, M, V, LDV ) )
+               LASTC = ILASLC( LASTV, N, C, LDC )
+*
+*              W := C**T * V**T  =  (C1**T * V1**T + C2**T * V2**T) (stored in WORK)
+*
+*              W := C2**T
+*
+               DO 190 J = 1, K
+                  CALL SCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
+     $                 WORK( 1, J ), 1 )
+  190          CONTINUE
+*
+*              W := W * V2**T
+*
+               CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1**T * V1**T
+*
+                  CALL SGEMM( 'Transpose', 'Transpose',
+     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**T  or  W * T
+*
+               CALL STRMM( 'Right', 'Lower', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V**T * W**T
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - V1**T * W**T
+*
+                  CALL SGEMM( 'Transpose', 'Transpose',
+     $                 LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2
+*
+               CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W**T
+*
+               DO 210 J = 1, K
+                  DO 200 I = 1, LASTC
+                     C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J)
+  200             CONTINUE
+  210          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**T  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILASLC( K, N, V, LDV ) )
+               LASTC = ILASLR( M, LASTV, C, LDC )
+*
+*              W := C * V**T  =  (C1*V1**T + C2*V2**T)  (stored in WORK)
+*
+*              W := C2
+*
+               DO 220 J = 1, K
+                  CALL SCOPY( LASTC, C( 1, LASTV-K+J ), 1,
+     $                 WORK( 1, J ), 1 )
+  220          CONTINUE
+*
+*              W := W * V2**T
+*
+               CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1 * V1**T
+*
+                  CALL SGEMM( 'No transpose', 'Transpose',
+     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**T
+*
+               CALL STRMM( 'Right', 'Lower', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - W * V1
+*
+                  CALL SGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2
+*
+               CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 240 J = 1, K
+                  DO 230 I = 1, LASTC
+                     C( I, LASTV-K+J ) = C( I, LASTV-K+J )
+     $                    - WORK( I, J )
+  230             CONTINUE
+  240          CONTINUE
+*
+            END IF
+*
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of SLARFB
+*
+      END
diff --git a/lapack/slarfg.f b/lapack/slarfg.f
new file mode 100644
index 000000000..4f10ffcaf
--- /dev/null
+++ b/lapack/slarfg.f
@@ -0,0 +1,196 @@
+*> \brief \b SLARFG
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download SLARFG + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarfg.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarfg.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarfg.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU )
+* 
+*       .. Scalar Arguments ..
+*       INTEGER            INCX, N
+*       REAL               ALPHA, TAU
+*       ..
+*       .. Array Arguments ..
+*       REAL               X( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> SLARFG generates a real elementary reflector H of order n, such
+*> that
+*>
+*>       H * ( alpha ) = ( beta ),   H**T * H = I.
+*>           (   x   )   (   0  )
+*>
+*> where alpha and beta are scalars, and x is an (n-1)-element real
+*> vector. H is represented in the form
+*>
+*>       H = I - tau * ( 1 ) * ( 1 v**T ) ,
+*>                     ( v )
+*>
+*> where tau is a real scalar and v is a real (n-1)-element
+*> vector.
+*>
+*> If the elements of x are all zero, then tau = 0 and H is taken to be
+*> the unit matrix.
+*>
+*> Otherwise  1 <= tau <= 2.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the elementary reflector.
+*> \endverbatim
+*>
+*> \param[in,out] ALPHA
+*> \verbatim
+*>          ALPHA is REAL
+*>          On entry, the value alpha.
+*>          On exit, it is overwritten with the value beta.
+*> \endverbatim
+*>
+*> \param[in,out] X
+*> \verbatim
+*>          X is REAL array, dimension
+*>                         (1+(N-2)*abs(INCX))
+*>          On entry, the vector x.
+*>          On exit, it is overwritten with the vector v.
+*> \endverbatim
+*>
+*> \param[in] INCX
+*> \verbatim
+*>          INCX is INTEGER
+*>          The increment between elements of X. INCX > 0.
+*> \endverbatim
+*>
+*> \param[out] TAU
+*> \verbatim
+*>          TAU is REAL
+*>          The value tau.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup realOTHERauxiliary
+*
+*  =====================================================================
+      SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      INTEGER            INCX, N
+      REAL               ALPHA, TAU
+*     ..
+*     .. Array Arguments ..
+      REAL               X( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE, ZERO
+      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            J, KNT
+      REAL               BETA, RSAFMN, SAFMIN, XNORM
+*     ..
+*     .. External Functions ..
+      REAL               SLAMCH, SLAPY2, SNRM2
+      EXTERNAL           SLAMCH, SLAPY2, SNRM2
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, SIGN
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SSCAL
+*     ..
+*     .. Executable Statements ..
+*
+      IF( N.LE.1 ) THEN
+         TAU = ZERO
+         RETURN
+      END IF
+*
+      XNORM = SNRM2( N-1, X, INCX )
+*
+      IF( XNORM.EQ.ZERO ) THEN
+*
+*        H  =  I
+*
+         TAU = ZERO
+      ELSE
+*
+*        general case
+*
+         BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
+         SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' )
+         KNT = 0
+         IF( ABS( BETA ).LT.SAFMIN ) THEN
+*
+*           XNORM, BETA may be inaccurate; scale X and recompute them
+*
+            RSAFMN = ONE / SAFMIN
+   10       CONTINUE
+            KNT = KNT + 1
+            CALL SSCAL( N-1, RSAFMN, X, INCX )
+            BETA = BETA*RSAFMN
+            ALPHA = ALPHA*RSAFMN
+            IF( ABS( BETA ).LT.SAFMIN )
+     $         GO TO 10
+*
+*           New BETA is at most 1, at least SAFMIN
+*
+            XNORM = SNRM2( N-1, X, INCX )
+            BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
+         END IF
+         TAU = ( BETA-ALPHA ) / BETA
+         CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
+*
+*        If ALPHA is subnormal, it may lose relative accuracy
+*
+         DO 20 J = 1, KNT
+            BETA = BETA*SAFMIN
+ 20      CONTINUE
+         ALPHA = BETA
+      END IF
+*
+      RETURN
+*
+*     End of SLARFG
+*
+      END
diff --git a/lapack/slarft.f b/lapack/slarft.f
new file mode 100644
index 000000000..30b0668e4
--- /dev/null
+++ b/lapack/slarft.f
@@ -0,0 +1,326 @@
+*> \brief \b SLARFT
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download SLARFT + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarft.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarft.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarft.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          DIRECT, STOREV
+*       INTEGER            K, LDT, LDV, N
+*       ..
+*       .. Array Arguments ..
+*       REAL               T( LDT, * ), TAU( * ), V( LDV, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> SLARFT forms the triangular factor T of a real block reflector H
+*> of order n, which is defined as a product of k elementary reflectors.
+*>
+*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
+*>
+*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
+*>
+*> If STOREV = 'C', the vector which defines the elementary reflector
+*> H(i) is stored in the i-th column of the array V, and
+*>
+*>    H  =  I - V * T * V**T
+*>
+*> If STOREV = 'R', the vector which defines the elementary reflector
+*> H(i) is stored in the i-th row of the array V, and
+*>
+*>    H  =  I - V**T * T * V
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] DIRECT
+*> \verbatim
+*>          DIRECT is CHARACTER*1
+*>          Specifies the order in which the elementary reflectors are
+*>          multiplied to form the block reflector:
+*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
+*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
+*> \endverbatim
+*>
+*> \param[in] STOREV
+*> \verbatim
+*>          STOREV is CHARACTER*1
+*>          Specifies how the vectors which define the elementary
+*>          reflectors are stored (see also Further Details):
+*>          = 'C': columnwise
+*>          = 'R': rowwise
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the block reflector H. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] K
+*> \verbatim
+*>          K is INTEGER
+*>          The order of the triangular factor T (= the number of
+*>          elementary reflectors). K >= 1.
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is REAL array, dimension
+*>                               (LDV,K) if STOREV = 'C'
+*>                               (LDV,N) if STOREV = 'R'
+*>          The matrix V. See further details.
+*> \endverbatim
+*>
+*> \param[in] LDV
+*> \verbatim
+*>          LDV is INTEGER
+*>          The leading dimension of the array V.
+*>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
+*> \endverbatim
+*>
+*> \param[in] TAU
+*> \verbatim
+*>          TAU is REAL array, dimension (K)
+*>          TAU(i) must contain the scalar factor of the elementary
+*>          reflector H(i).
+*> \endverbatim
+*>
+*> \param[out] T
+*> \verbatim
+*>          T is REAL array, dimension (LDT,K)
+*>          The k by k triangular factor T of the block reflector.
+*>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
+*>          lower triangular. The rest of the array is not used.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*>          LDT is INTEGER
+*>          The leading dimension of the array T. LDT >= K.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date April 2012
+*
+*> \ingroup realOTHERauxiliary
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  The shape of the matrix V and the storage of the vectors which define
+*>  the H(i) is best illustrated by the following example with n = 5 and
+*>  k = 3. The elements equal to 1 are not stored.
+*>
+*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
+*>
+*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
+*>                   ( v1  1    )                     (     1 v2 v2 v2 )
+*>                   ( v1 v2  1 )                     (        1 v3 v3 )
+*>                   ( v1 v2 v3 )
+*>                   ( v1 v2 v3 )
+*>
+*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
+*>
+*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
+*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
+*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
+*>                   (     1 v3 )
+*>                   (        1 )
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
+*
+*  -- LAPACK auxiliary routine (version 3.4.1) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     April 2012
+*
+*     .. Scalar Arguments ..
+      CHARACTER          DIRECT, STOREV
+      INTEGER            K, LDT, LDV, N
+*     ..
+*     .. Array Arguments ..
+      REAL               T( LDT, * ), TAU( * ), V( LDV, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      REAL               ONE, ZERO
+      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, J, PREVLASTV, LASTV
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           SGEMV, STRMV
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      IF( LSAME( DIRECT, 'F' ) ) THEN
+         PREVLASTV = N
+         DO I = 1, K
+            PREVLASTV = MAX( I, PREVLASTV )
+            IF( TAU( I ).EQ.ZERO ) THEN
+*
+*              H(i)  =  I
+*
+               DO J = 1, I
+                  T( J, I ) = ZERO
+               END DO
+            ELSE
+*
+*              general case
+*
+               IF( LSAME( STOREV, 'C' ) ) THEN
+*                 Skip any trailing zeros.
+                  DO LASTV = N, I+1, -1
+                     IF( V( LASTV, I ).NE.ZERO ) EXIT
+                  END DO
+                  DO J = 1, I-1
+                     T( J, I ) = -TAU( I ) * V( I , J )
+                  END DO   
+                  J = MIN( LASTV, PREVLASTV )
+*
+*                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i)
+*
+                  CALL SGEMV( 'Transpose', J-I, I-1, -TAU( I ),
+     $                        V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE,
+     $                        T( 1, I ), 1 )
+               ELSE
+*                 Skip any trailing zeros.
+                  DO LASTV = N, I+1, -1
+                     IF( V( I, LASTV ).NE.ZERO ) EXIT
+                  END DO
+                  DO J = 1, I-1
+                     T( J, I ) = -TAU( I ) * V( J , I )
+                  END DO   
+                  J = MIN( LASTV, PREVLASTV )
+*
+*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T
+*
+                  CALL SGEMV( 'No transpose', I-1, J-I, -TAU( I ),
+     $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV, 
+     $                        ONE, T( 1, I ), 1 )
+               END IF
+*
+*              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
+*
+               CALL STRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
+     $                     LDT, T( 1, I ), 1 )
+               T( I, I ) = TAU( I )
+               IF( I.GT.1 ) THEN
+                  PREVLASTV = MAX( PREVLASTV, LASTV )
+               ELSE
+                  PREVLASTV = LASTV
+               END IF
+            END IF
+         END DO
+      ELSE
+         PREVLASTV = 1
+         DO I = K, 1, -1
+            IF( TAU( I ).EQ.ZERO ) THEN
+*
+*              H(i)  =  I
+*
+               DO J = I, K
+                  T( J, I ) = ZERO
+               END DO
+            ELSE
+*
+*              general case
+*
+               IF( I.LT.K ) THEN
+                  IF( LSAME( STOREV, 'C' ) ) THEN
+*                    Skip any leading zeros.
+                     DO LASTV = 1, I-1
+                        IF( V( LASTV, I ).NE.ZERO ) EXIT
+                     END DO
+                     DO J = I+1, K
+                        T( J, I ) = -TAU( I ) * V( N-K+I , J )
+                     END DO   
+                     J = MAX( LASTV, PREVLASTV )
+*
+*                    T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i)
+*
+                     CALL SGEMV( 'Transpose', N-K+I-J, K-I, -TAU( I ),
+     $                           V( J, I+1 ), LDV, V( J, I ), 1, ONE,
+     $                           T( I+1, I ), 1 )
+                  ELSE
+*                    Skip any leading zeros.
+                     DO LASTV = 1, I-1
+                        IF( V( I, LASTV ).NE.ZERO ) EXIT
+                     END DO
+                     DO J = I+1, K
+                        T( J, I ) = -TAU( I ) * V( J, N-K+I )
+                     END DO   
+                     J = MAX( LASTV, PREVLASTV )
+*
+*                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T
+*
+                     CALL SGEMV( 'No transpose', K-I, N-K+I-J,
+     $                    -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV,
+     $                    ONE, T( I+1, I ), 1 )
+                  END IF
+*
+*                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
+*
+                  CALL STRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
+     $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
+                  IF( I.GT.1 ) THEN
+                     PREVLASTV = MIN( PREVLASTV, LASTV )
+                  ELSE
+                     PREVLASTV = LASTV
+                  END IF
+               END IF
+               T( I, I ) = TAU( I )
+            END IF
+         END DO
+      END IF
+      RETURN
+*
+*     End of SLARFT
+*
+      END
diff --git a/lapack/zlarf.f b/lapack/zlarf.f
new file mode 100644
index 000000000..53f314d64
--- /dev/null
+++ b/lapack/zlarf.f
@@ -0,0 +1,232 @@
+*> \brief \b ZLARF
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download ZLARF + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarf.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarf.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarf.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          SIDE
+*       INTEGER            INCV, LDC, M, N
+*       COMPLEX*16         TAU
+*       ..
+*       .. Array Arguments ..
+*       COMPLEX*16         C( LDC, * ), V( * ), WORK( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> ZLARF applies a complex elementary reflector H to a complex M-by-N
+*> matrix C, from either the left or the right. H is represented in the
+*> form
+*>
+*>       H = I - tau * v * v**H
+*>
+*> where tau is a complex scalar and v is a complex vector.
+*>
+*> If tau = 0, then H is taken to be the unit matrix.
+*>
+*> To apply H**H, supply conjg(tau) instead
+*> tau.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*>          SIDE is CHARACTER*1
+*>          = 'L': form  H * C
+*>          = 'R': form  C * H
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*>          M is INTEGER
+*>          The number of rows of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of columns of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is COMPLEX*16 array, dimension
+*>                     (1 + (M-1)*abs(INCV)) if SIDE = 'L'
+*>                  or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
+*>          The vector v in the representation of H. V is not used if
+*>          TAU = 0.
+*> \endverbatim
+*>
+*> \param[in] INCV
+*> \verbatim
+*>          INCV is INTEGER
+*>          The increment between elements of v. INCV <> 0.
+*> \endverbatim
+*>
+*> \param[in] TAU
+*> \verbatim
+*>          TAU is COMPLEX*16
+*>          The value tau in the representation of H.
+*> \endverbatim
+*>
+*> \param[in,out] C
+*> \verbatim
+*>          C is COMPLEX*16 array, dimension (LDC,N)
+*>          On entry, the M-by-N matrix C.
+*>          On exit, C is overwritten by the matrix H * C if SIDE = 'L',
+*>          or C * H if SIDE = 'R'.
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX*16 array, dimension
+*>                         (N) if SIDE = 'L'
+*>                      or (M) if SIDE = 'R'
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complex16OTHERauxiliary
+*
+*  =====================================================================
+      SUBROUTINE ZLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          SIDE
+      INTEGER            INCV, LDC, M, N
+      COMPLEX*16         TAU
+*     ..
+*     .. Array Arguments ..
+      COMPLEX*16         C( LDC, * ), V( * ), WORK( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX*16         ONE, ZERO
+      PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
+     $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
+*     ..
+*     .. Local Scalars ..
+      LOGICAL            APPLYLEFT
+      INTEGER            I, LASTV, LASTC
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZGEMV, ZGERC
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAZLR, ILAZLC
+      EXTERNAL           LSAME, ILAZLR, ILAZLC
+*     ..
+*     .. Executable Statements ..
+*
+      APPLYLEFT = LSAME( SIDE, 'L' )
+      LASTV = 0
+      LASTC = 0
+      IF( TAU.NE.ZERO ) THEN
+*     Set up variables for scanning V.  LASTV begins pointing to the end
+*     of V.
+         IF( APPLYLEFT ) THEN
+            LASTV = M
+         ELSE
+            LASTV = N
+         END IF
+         IF( INCV.GT.0 ) THEN
+            I = 1 + (LASTV-1) * INCV
+         ELSE
+            I = 1
+         END IF
+*     Look for the last non-zero row in V.
+         DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO )
+            LASTV = LASTV - 1
+            I = I - INCV
+         END DO
+         IF( APPLYLEFT ) THEN
+*     Scan for the last non-zero column in C(1:lastv,:).
+            LASTC = ILAZLC(LASTV, N, C, LDC)
+         ELSE
+*     Scan for the last non-zero row in C(:,1:lastv).
+            LASTC = ILAZLR(M, LASTV, C, LDC)
+         END IF
+      END IF
+*     Note that lastc.eq.0 renders the BLAS operations null; no special
+*     case is needed at this level.
+      IF( APPLYLEFT ) THEN
+*
+*        Form  H * C
+*
+         IF( LASTV.GT.0 ) THEN
+*
+*           w(1:lastc,1) := C(1:lastv,1:lastc)**H * v(1:lastv,1)
+*
+            CALL ZGEMV( 'Conjugate transpose', LASTV, LASTC, ONE,
+     $           C, LDC, V, INCV, ZERO, WORK, 1 )
+*
+*           C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)**H
+*
+            CALL ZGERC( LASTV, LASTC, -TAU, V, INCV, WORK, 1, C, LDC )
+         END IF
+      ELSE
+*
+*        Form  C * H
+*
+         IF( LASTV.GT.0 ) THEN
+*
+*           w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1)
+*
+            CALL ZGEMV( 'No transpose', LASTC, LASTV, ONE, C, LDC,
+     $           V, INCV, ZERO, WORK, 1 )
+*
+*           C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)**H
+*
+            CALL ZGERC( LASTC, LASTV, -TAU, WORK, 1, V, INCV, C, LDC )
+         END IF
+      END IF
+      RETURN
+*
+*     End of ZLARF
+*
+      END
diff --git a/lapack/zlarfb.f b/lapack/zlarfb.f
new file mode 100644
index 000000000..30fc4b940
--- /dev/null
+++ b/lapack/zlarfb.f
@@ -0,0 +1,774 @@
+*> \brief \b ZLARFB
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download ZLARFB + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarfb.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarfb.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarfb.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
+*                          T, LDT, C, LDC, WORK, LDWORK )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          DIRECT, SIDE, STOREV, TRANS
+*       INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
+*       ..
+*       .. Array Arguments ..
+*       COMPLEX*16         C( LDC, * ), T( LDT, * ), V( LDV, * ),
+*      $                   WORK( LDWORK, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> ZLARFB applies a complex block reflector H or its transpose H**H to a
+*> complex M-by-N matrix C, from either the left or the right.
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] SIDE
+*> \verbatim
+*>          SIDE is CHARACTER*1
+*>          = 'L': apply H or H**H from the Left
+*>          = 'R': apply H or H**H from the Right
+*> \endverbatim
+*>
+*> \param[in] TRANS
+*> \verbatim
+*>          TRANS is CHARACTER*1
+*>          = 'N': apply H (No transpose)
+*>          = 'C': apply H**H (Conjugate transpose)
+*> \endverbatim
+*>
+*> \param[in] DIRECT
+*> \verbatim
+*>          DIRECT is CHARACTER*1
+*>          Indicates how H is formed from a product of elementary
+*>          reflectors
+*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
+*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
+*> \endverbatim
+*>
+*> \param[in] STOREV
+*> \verbatim
+*>          STOREV is CHARACTER*1
+*>          Indicates how the vectors which define the elementary
+*>          reflectors are stored:
+*>          = 'C': Columnwise
+*>          = 'R': Rowwise
+*> \endverbatim
+*>
+*> \param[in] M
+*> \verbatim
+*>          M is INTEGER
+*>          The number of rows of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The number of columns of the matrix C.
+*> \endverbatim
+*>
+*> \param[in] K
+*> \verbatim
+*>          K is INTEGER
+*>          The order of the matrix T (= the number of elementary
+*>          reflectors whose product defines the block reflector).
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is COMPLEX*16 array, dimension
+*>                                (LDV,K) if STOREV = 'C'
+*>                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
+*>                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
+*>          See Further Details.
+*> \endverbatim
+*>
+*> \param[in] LDV
+*> \verbatim
+*>          LDV is INTEGER
+*>          The leading dimension of the array V.
+*>          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
+*>          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
+*>          if STOREV = 'R', LDV >= K.
+*> \endverbatim
+*>
+*> \param[in] T
+*> \verbatim
+*>          T is COMPLEX*16 array, dimension (LDT,K)
+*>          The triangular K-by-K matrix T in the representation of the
+*>          block reflector.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*>          LDT is INTEGER
+*>          The leading dimension of the array T. LDT >= K.
+*> \endverbatim
+*>
+*> \param[in,out] C
+*> \verbatim
+*>          C is COMPLEX*16 array, dimension (LDC,N)
+*>          On entry, the M-by-N matrix C.
+*>          On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H.
+*> \endverbatim
+*>
+*> \param[in] LDC
+*> \verbatim
+*>          LDC is INTEGER
+*>          The leading dimension of the array C. LDC >= max(1,M).
+*> \endverbatim
+*>
+*> \param[out] WORK
+*> \verbatim
+*>          WORK is COMPLEX*16 array, dimension (LDWORK,K)
+*> \endverbatim
+*>
+*> \param[in] LDWORK
+*> \verbatim
+*>          LDWORK is INTEGER
+*>          The leading dimension of the array WORK.
+*>          If SIDE = 'L', LDWORK >= max(1,N);
+*>          if SIDE = 'R', LDWORK >= max(1,M).
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complex16OTHERauxiliary
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  The shape of the matrix V and the storage of the vectors which define
+*>  the H(i) is best illustrated by the following example with n = 5 and
+*>  k = 3. The elements equal to 1 are not stored; the corresponding
+*>  array elements are modified but restored on exit. The rest of the
+*>  array is not used.
+*>
+*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
+*>
+*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
+*>                   ( v1  1    )                     (     1 v2 v2 v2 )
+*>                   ( v1 v2  1 )                     (        1 v3 v3 )
+*>                   ( v1 v2 v3 )
+*>                   ( v1 v2 v3 )
+*>
+*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
+*>
+*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
+*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
+*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
+*>                   (     1 v3 )
+*>                   (        1 )
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE ZLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV,
+     $                   T, LDT, C, LDC, WORK, LDWORK )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      CHARACTER          DIRECT, SIDE, STOREV, TRANS
+      INTEGER            K, LDC, LDT, LDV, LDWORK, M, N
+*     ..
+*     .. Array Arguments ..
+      COMPLEX*16         C( LDC, * ), T( LDT, * ), V( LDV, * ),
+     $                   WORK( LDWORK, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX*16         ONE
+      PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
+*     ..
+*     .. Local Scalars ..
+      CHARACTER          TRANST
+      INTEGER            I, J, LASTV, LASTC
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      INTEGER            ILAZLR, ILAZLC
+      EXTERNAL           LSAME, ILAZLR, ILAZLC
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZCOPY, ZGEMM, ZLACGV, ZTRMM
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          DCONJG
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick return if possible
+*
+      IF( M.LE.0 .OR. N.LE.0 )
+     $   RETURN
+*
+      IF( LSAME( TRANS, 'N' ) ) THEN
+         TRANST = 'C'
+      ELSE
+         TRANST = 'N'
+      END IF
+*
+      IF( LSAME( STOREV, 'C' ) ) THEN
+*
+         IF( LSAME( DIRECT, 'F' ) ) THEN
+*
+*           Let  V =  ( V1 )    (first K rows)
+*                     ( V2 )
+*           where  V1  is unit lower triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**H * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILAZLR( M, K, V, LDV ) )
+               LASTC = ILAZLC( LASTV, N, C, LDC )
+*
+*              W := C**H * V  =  (C1**H * V1 + C2**H * V2)  (stored in WORK)
+*
+*              W := C1**H
+*
+               DO 10 J = 1, K
+                  CALL ZCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
+                  CALL ZLACGV( LASTC, WORK( 1, J ), 1 )
+   10          CONTINUE
+*
+*              W := W * V1
+*
+               CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2**H *V2
+*
+                  CALL ZGEMM( 'Conjugate transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K, ONE, C( K+1, 1 ), LDC,
+     $                 V( K+1, 1 ), LDV, ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**H  or  W * T
+*
+               CALL ZTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V * W**H
+*
+               IF( M.GT.K ) THEN
+*
+*                 C2 := C2 - V2 * W**H
+*
+                  CALL ZGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTV-K, LASTC, K,
+     $                 -ONE, V( K+1, 1 ), LDV, WORK, LDWORK,
+     $                 ONE, C( K+1, 1 ), LDC )
+               END IF
+*
+*              W := W * V1**H
+*
+               CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W**H
+*
+               DO 30 J = 1, K
+                  DO 20 I = 1, LASTC
+                     C( J, I ) = C( J, I ) - DCONJG( WORK( I, J ) )
+   20             CONTINUE
+   30          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILAZLR( N, K, V, LDV ) )
+               LASTC = ILAZLR( M, LASTV, C, LDC )
+*
+*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
+*
+*              W := C1
+*
+               DO 40 J = 1, K
+                  CALL ZCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
+   40          CONTINUE
+*
+*              W := W * V1
+*
+               CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2 * V2
+*
+                  CALL ZGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**H
+*
+               CALL ZTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V**H
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - W * V2**H
+*
+                  CALL ZGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTC, LASTV-K, K,
+     $                 -ONE, WORK, LDWORK, V( K+1, 1 ), LDV,
+     $                 ONE, C( 1, K+1 ), LDC )
+               END IF
+*
+*              W := W * V1**H
+*
+               CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 60 J = 1, K
+                  DO 50 I = 1, LASTC
+                     C( I, J ) = C( I, J ) - WORK( I, J )
+   50             CONTINUE
+   60          CONTINUE
+            END IF
+*
+         ELSE
+*
+*           Let  V =  ( V1 )
+*                     ( V2 )    (last K rows)
+*           where  V2  is unit upper triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**H * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILAZLR( M, K, V, LDV ) )
+               LASTC = ILAZLC( LASTV, N, C, LDC )
+*
+*              W := C**H * V  =  (C1**H * V1 + C2**H * V2)  (stored in WORK)
+*
+*              W := C2**H
+*
+               DO 70 J = 1, K
+                  CALL ZCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
+     $                 WORK( 1, J ), 1 )
+                  CALL ZLACGV( LASTC, WORK( 1, J ), 1 )
+   70          CONTINUE
+*
+*              W := W * V2
+*
+               CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1**H*V1
+*
+                  CALL ZGEMM( 'Conjugate transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C, LDC, V, LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**H  or  W * T
+*
+               CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V * W**H
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - V1 * W**H
+*
+                  CALL ZGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTV-K, LASTC, K,
+     $                 -ONE, V, LDV, WORK, LDWORK,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2**H
+*
+               CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W**H
+*
+               DO 90 J = 1, K
+                  DO 80 I = 1, LASTC
+                     C( LASTV-K+J, I ) = C( LASTV-K+J, I ) -
+     $                               DCONJG( WORK( I, J ) )
+   80             CONTINUE
+   90          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILAZLR( N, K, V, LDV ) )
+               LASTC = ILAZLR( M, LASTV, C, LDC )
+*
+*              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
+*
+*              W := C2
+*
+               DO 100 J = 1, K
+                  CALL ZCOPY( LASTC, C( 1, LASTV-K+J ), 1,
+     $                 WORK( 1, J ), 1 )
+  100          CONTINUE
+*
+*              W := W * V2
+*
+               CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1 * V1
+*
+                  CALL ZGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, K, LASTV-K,
+     $                 ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**H
+*
+               CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V**H
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - W * V1**H
+*
+                  CALL ZGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2**H
+*
+               CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W
+*
+               DO 120 J = 1, K
+                  DO 110 I = 1, LASTC
+                     C( I, LASTV-K+J ) = C( I, LASTV-K+J )
+     $                    - WORK( I, J )
+  110             CONTINUE
+  120          CONTINUE
+            END IF
+         END IF
+*
+      ELSE IF( LSAME( STOREV, 'R' ) ) THEN
+*
+         IF( LSAME( DIRECT, 'F' ) ) THEN
+*
+*           Let  V =  ( V1  V2 )    (V1: first K columns)
+*           where  V1  is unit upper triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**H * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILAZLC( K, M, V, LDV ) )
+               LASTC = ILAZLC( LASTV, N, C, LDC )
+*
+*              W := C**H * V**H  =  (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
+*
+*              W := C1**H
+*
+               DO 130 J = 1, K
+                  CALL ZCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 )
+                  CALL ZLACGV( LASTC, WORK( 1, J ), 1 )
+  130          CONTINUE
+*
+*              W := W * V1**H
+*
+               CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
+     $                     'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2**H*V2**H
+*
+                  CALL ZGEMM( 'Conjugate transpose',
+     $                 'Conjugate transpose', LASTC, K, LASTV-K,
+     $                 ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV,
+     $                 ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**H  or  W * T
+*
+               CALL ZTRMM( 'Right', 'Upper', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V**H * W**H
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - V2**H * W**H
+*
+                  CALL ZGEMM( 'Conjugate transpose',
+     $                 'Conjugate transpose', LASTV-K, LASTC, K,
+     $                 -ONE, V( 1, K+1 ), LDV, WORK, LDWORK,
+     $                 ONE, C( K+1, 1 ), LDC )
+               END IF
+*
+*              W := W * V1
+*
+               CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W**H
+*
+               DO 150 J = 1, K
+                  DO 140 I = 1, LASTC
+                     C( J, I ) = C( J, I ) - DCONJG( WORK( I, J ) )
+  140             CONTINUE
+  150          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILAZLC( K, N, V, LDV ) )
+               LASTC = ILAZLR( M, LASTV, C, LDC )
+*
+*              W := C * V**H  =  (C1*V1**H + C2*V2**H)  (stored in WORK)
+*
+*              W := C1
+*
+               DO 160 J = 1, K
+                  CALL ZCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 )
+  160          CONTINUE
+*
+*              W := W * V1**H
+*
+               CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose',
+     $                     'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C2 * V2**H
+*
+                  CALL ZGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTC, K, LASTV-K, ONE, C( 1, K+1 ), LDC,
+     $                 V( 1, K+1 ), LDV, ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**H
+*
+               CALL ZTRMM( 'Right', 'Upper', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C2 := C2 - W * V2
+*
+                  CALL ZGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, LASTV-K, K,
+     $                 -ONE, WORK, LDWORK, V( 1, K+1 ), LDV,
+     $                 ONE, C( 1, K+1 ), LDC )
+               END IF
+*
+*              W := W * V1
+*
+               CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V, LDV, WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 180 J = 1, K
+                  DO 170 I = 1, LASTC
+                     C( I, J ) = C( I, J ) - WORK( I, J )
+  170             CONTINUE
+  180          CONTINUE
+*
+            END IF
+*
+         ELSE
+*
+*           Let  V =  ( V1  V2 )    (V2: last K columns)
+*           where  V2  is unit lower triangular.
+*
+            IF( LSAME( SIDE, 'L' ) ) THEN
+*
+*              Form  H * C  or  H**H * C  where  C = ( C1 )
+*                                                    ( C2 )
+*
+               LASTV = MAX( K, ILAZLC( K, M, V, LDV ) )
+               LASTC = ILAZLC( LASTV, N, C, LDC )
+*
+*              W := C**H * V**H  =  (C1**H * V1**H + C2**H * V2**H) (stored in WORK)
+*
+*              W := C2**H
+*
+               DO 190 J = 1, K
+                  CALL ZCOPY( LASTC, C( LASTV-K+J, 1 ), LDC,
+     $                 WORK( 1, J ), 1 )
+                  CALL ZLACGV( LASTC, WORK( 1, J ), 1 )
+  190          CONTINUE
+*
+*              W := W * V2**H
+*
+               CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1**H * V1**H
+*
+                  CALL ZGEMM( 'Conjugate transpose',
+     $                 'Conjugate transpose', LASTC, K, LASTV-K,
+     $                 ONE, C, LDC, V, LDV, ONE, WORK, LDWORK )
+               END IF
+*
+*              W := W * T**H  or  W * T
+*
+               CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - V**H * W**H
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - V1**H * W**H
+*
+                  CALL ZGEMM( 'Conjugate transpose',
+     $                 'Conjugate transpose', LASTV-K, LASTC, K,
+     $                 -ONE, V, LDV, WORK, LDWORK, ONE, C, LDC )
+               END IF
+*
+*              W := W * V2
+*
+               CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C2 := C2 - W**H
+*
+               DO 210 J = 1, K
+                  DO 200 I = 1, LASTC
+                     C( LASTV-K+J, I ) = C( LASTV-K+J, I ) -
+     $                               DCONJG( WORK( I, J ) )
+  200             CONTINUE
+  210          CONTINUE
+*
+            ELSE IF( LSAME( SIDE, 'R' ) ) THEN
+*
+*              Form  C * H  or  C * H**H  where  C = ( C1  C2 )
+*
+               LASTV = MAX( K, ILAZLC( K, N, V, LDV ) )
+               LASTC = ILAZLR( M, LASTV, C, LDC )
+*
+*              W := C * V**H  =  (C1*V1**H + C2*V2**H)  (stored in WORK)
+*
+*              W := C2
+*
+               DO 220 J = 1, K
+                  CALL ZCOPY( LASTC, C( 1, LASTV-K+J ), 1,
+     $                 WORK( 1, J ), 1 )
+  220          CONTINUE
+*
+*              W := W * V2**H
+*
+               CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
+     $              'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+               IF( LASTV.GT.K ) THEN
+*
+*                 W := W + C1 * V1**H
+*
+                  CALL ZGEMM( 'No transpose', 'Conjugate transpose',
+     $                 LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, ONE,
+     $                 WORK, LDWORK )
+               END IF
+*
+*              W := W * T  or  W * T**H
+*
+               CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit',
+     $              LASTC, K, ONE, T, LDT, WORK, LDWORK )
+*
+*              C := C - W * V
+*
+               IF( LASTV.GT.K ) THEN
+*
+*                 C1 := C1 - W * V1
+*
+                  CALL ZGEMM( 'No transpose', 'No transpose',
+     $                 LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV,
+     $                 ONE, C, LDC )
+               END IF
+*
+*              W := W * V2
+*
+               CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit',
+     $              LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV,
+     $              WORK, LDWORK )
+*
+*              C1 := C1 - W
+*
+               DO 240 J = 1, K
+                  DO 230 I = 1, LASTC
+                     C( I, LASTV-K+J ) = C( I, LASTV-K+J )
+     $                    - WORK( I, J )
+  230             CONTINUE
+  240          CONTINUE
+*
+            END IF
+*
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of ZLARFB
+*
+      END
diff --git a/lapack/zlarfg.f b/lapack/zlarfg.f
new file mode 100644
index 000000000..a90ae9f74
--- /dev/null
+++ b/lapack/zlarfg.f
@@ -0,0 +1,203 @@
+*> \brief \b ZLARFG
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download ZLARFG + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarfg.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarfg.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarfg.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU )
+* 
+*       .. Scalar Arguments ..
+*       INTEGER            INCX, N
+*       COMPLEX*16         ALPHA, TAU
+*       ..
+*       .. Array Arguments ..
+*       COMPLEX*16         X( * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> ZLARFG generates a complex elementary reflector H of order n, such
+*> that
+*>
+*>       H**H * ( alpha ) = ( beta ),   H**H * H = I.
+*>              (   x   )   (   0  )
+*>
+*> where alpha and beta are scalars, with beta real, and x is an
+*> (n-1)-element complex vector. H is represented in the form
+*>
+*>       H = I - tau * ( 1 ) * ( 1 v**H ) ,
+*>                     ( v )
+*>
+*> where tau is a complex scalar and v is a complex (n-1)-element
+*> vector. Note that H is not hermitian.
+*>
+*> If the elements of x are all zero and alpha is real, then tau = 0
+*> and H is taken to be the unit matrix.
+*>
+*> Otherwise  1 <= real(tau) <= 2  and  abs(tau-1) <= 1 .
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the elementary reflector.
+*> \endverbatim
+*>
+*> \param[in,out] ALPHA
+*> \verbatim
+*>          ALPHA is COMPLEX*16
+*>          On entry, the value alpha.
+*>          On exit, it is overwritten with the value beta.
+*> \endverbatim
+*>
+*> \param[in,out] X
+*> \verbatim
+*>          X is COMPLEX*16 array, dimension
+*>                         (1+(N-2)*abs(INCX))
+*>          On entry, the vector x.
+*>          On exit, it is overwritten with the vector v.
+*> \endverbatim
+*>
+*> \param[in] INCX
+*> \verbatim
+*>          INCX is INTEGER
+*>          The increment between elements of X. INCX > 0.
+*> \endverbatim
+*>
+*> \param[out] TAU
+*> \verbatim
+*>          TAU is COMPLEX*16
+*>          The value tau.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date November 2011
+*
+*> \ingroup complex16OTHERauxiliary
+*
+*  =====================================================================
+      SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU )
+*
+*  -- LAPACK auxiliary routine (version 3.4.0) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     November 2011
+*
+*     .. Scalar Arguments ..
+      INTEGER            INCX, N
+      COMPLEX*16         ALPHA, TAU
+*     ..
+*     .. Array Arguments ..
+      COMPLEX*16         X( * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      DOUBLE PRECISION   ONE, ZERO
+      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            J, KNT
+      DOUBLE PRECISION   ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM
+*     ..
+*     .. External Functions ..
+      DOUBLE PRECISION   DLAMCH, DLAPY3, DZNRM2
+      COMPLEX*16         ZLADIV
+      EXTERNAL           DLAMCH, DLAPY3, DZNRM2, ZLADIV
+*     ..
+*     .. Intrinsic Functions ..
+      INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, SIGN
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZDSCAL, ZSCAL
+*     ..
+*     .. Executable Statements ..
+*
+      IF( N.LE.0 ) THEN
+         TAU = ZERO
+         RETURN
+      END IF
+*
+      XNORM = DZNRM2( N-1, X, INCX )
+      ALPHR = DBLE( ALPHA )
+      ALPHI = DIMAG( ALPHA )
+*
+      IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN
+*
+*        H  =  I
+*
+         TAU = ZERO
+      ELSE
+*
+*        general case
+*
+         BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
+         SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' )
+         RSAFMN = ONE / SAFMIN
+*
+         KNT = 0
+         IF( ABS( BETA ).LT.SAFMIN ) THEN
+*
+*           XNORM, BETA may be inaccurate; scale X and recompute them
+*
+   10       CONTINUE
+            KNT = KNT + 1
+            CALL ZDSCAL( N-1, RSAFMN, X, INCX )
+            BETA = BETA*RSAFMN
+            ALPHI = ALPHI*RSAFMN
+            ALPHR = ALPHR*RSAFMN
+            IF( ABS( BETA ).LT.SAFMIN )
+     $         GO TO 10
+*
+*           New BETA is at most 1, at least SAFMIN
+*
+            XNORM = DZNRM2( N-1, X, INCX )
+            ALPHA = DCMPLX( ALPHR, ALPHI )
+            BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
+         END IF
+         TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
+         ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA )
+         CALL ZSCAL( N-1, ALPHA, X, INCX )
+*
+*        If ALPHA is subnormal, it may lose relative accuracy
+*
+         DO 20 J = 1, KNT
+            BETA = BETA*SAFMIN
+ 20      CONTINUE
+         ALPHA = BETA
+      END IF
+*
+      RETURN
+*
+*     End of ZLARFG
+*
+      END
diff --git a/lapack/zlarft.f b/lapack/zlarft.f
new file mode 100644
index 000000000..6a6151fd0
--- /dev/null
+++ b/lapack/zlarft.f
@@ -0,0 +1,327 @@
+*> \brief \b ZLARFT
+*
+*  =========== DOCUMENTATION ===========
+*
+* Online html documentation available at 
+*            http://www.netlib.org/lapack/explore-html/ 
+*
+*> \htmlonly
+*> Download ZLARFT + dependencies 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarft.f"> 
+*> [TGZ]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarft.f"> 
+*> [ZIP]</a> 
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarft.f"> 
+*> [TXT]</a>
+*> \endhtmlonly 
+*
+*  Definition:
+*  ===========
+*
+*       SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
+* 
+*       .. Scalar Arguments ..
+*       CHARACTER          DIRECT, STOREV
+*       INTEGER            K, LDT, LDV, N
+*       ..
+*       .. Array Arguments ..
+*       COMPLEX*16         T( LDT, * ), TAU( * ), V( LDV, * )
+*       ..
+*  
+*
+*> \par Purpose:
+*  =============
+*>
+*> \verbatim
+*>
+*> ZLARFT forms the triangular factor T of a complex block reflector H
+*> of order n, which is defined as a product of k elementary reflectors.
+*>
+*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
+*>
+*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
+*>
+*> If STOREV = 'C', the vector which defines the elementary reflector
+*> H(i) is stored in the i-th column of the array V, and
+*>
+*>    H  =  I - V * T * V**H
+*>
+*> If STOREV = 'R', the vector which defines the elementary reflector
+*> H(i) is stored in the i-th row of the array V, and
+*>
+*>    H  =  I - V**H * T * V
+*> \endverbatim
+*
+*  Arguments:
+*  ==========
+*
+*> \param[in] DIRECT
+*> \verbatim
+*>          DIRECT is CHARACTER*1
+*>          Specifies the order in which the elementary reflectors are
+*>          multiplied to form the block reflector:
+*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
+*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
+*> \endverbatim
+*>
+*> \param[in] STOREV
+*> \verbatim
+*>          STOREV is CHARACTER*1
+*>          Specifies how the vectors which define the elementary
+*>          reflectors are stored (see also Further Details):
+*>          = 'C': columnwise
+*>          = 'R': rowwise
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*>          N is INTEGER
+*>          The order of the block reflector H. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] K
+*> \verbatim
+*>          K is INTEGER
+*>          The order of the triangular factor T (= the number of
+*>          elementary reflectors). K >= 1.
+*> \endverbatim
+*>
+*> \param[in] V
+*> \verbatim
+*>          V is COMPLEX*16 array, dimension
+*>                               (LDV,K) if STOREV = 'C'
+*>                               (LDV,N) if STOREV = 'R'
+*>          The matrix V. See further details.
+*> \endverbatim
+*>
+*> \param[in] LDV
+*> \verbatim
+*>          LDV is INTEGER
+*>          The leading dimension of the array V.
+*>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
+*> \endverbatim
+*>
+*> \param[in] TAU
+*> \verbatim
+*>          TAU is COMPLEX*16 array, dimension (K)
+*>          TAU(i) must contain the scalar factor of the elementary
+*>          reflector H(i).
+*> \endverbatim
+*>
+*> \param[out] T
+*> \verbatim
+*>          T is COMPLEX*16 array, dimension (LDT,K)
+*>          The k by k triangular factor T of the block reflector.
+*>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
+*>          lower triangular. The rest of the array is not used.
+*> \endverbatim
+*>
+*> \param[in] LDT
+*> \verbatim
+*>          LDT is INTEGER
+*>          The leading dimension of the array T. LDT >= K.
+*> \endverbatim
+*
+*  Authors:
+*  ========
+*
+*> \author Univ. of Tennessee 
+*> \author Univ. of California Berkeley 
+*> \author Univ. of Colorado Denver 
+*> \author NAG Ltd. 
+*
+*> \date April 2012
+*
+*> \ingroup complex16OTHERauxiliary
+*
+*> \par Further Details:
+*  =====================
+*>
+*> \verbatim
+*>
+*>  The shape of the matrix V and the storage of the vectors which define
+*>  the H(i) is best illustrated by the following example with n = 5 and
+*>  k = 3. The elements equal to 1 are not stored.
+*>
+*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
+*>
+*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
+*>                   ( v1  1    )                     (     1 v2 v2 v2 )
+*>                   ( v1 v2  1 )                     (        1 v3 v3 )
+*>                   ( v1 v2 v3 )
+*>                   ( v1 v2 v3 )
+*>
+*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
+*>
+*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
+*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
+*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
+*>                   (     1 v3 )
+*>                   (        1 )
+*> \endverbatim
+*>
+*  =====================================================================
+      SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
+*
+*  -- LAPACK auxiliary routine (version 3.4.1) --
+*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
+*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+*     April 2012
+*
+*     .. Scalar Arguments ..
+      CHARACTER          DIRECT, STOREV
+      INTEGER            K, LDT, LDV, N
+*     ..
+*     .. Array Arguments ..
+      COMPLEX*16         T( LDT, * ), TAU( * ), V( LDV, * )
+*     ..
+*
+*  =====================================================================
+*
+*     .. Parameters ..
+      COMPLEX*16         ONE, ZERO
+      PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
+     $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
+*     ..
+*     .. Local Scalars ..
+      INTEGER            I, J, PREVLASTV, LASTV
+*     ..
+*     .. External Subroutines ..
+      EXTERNAL           ZGEMV, ZLACGV, ZTRMV
+*     ..
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     ..
+*     .. Executable Statements ..
+*
+*     Quick return if possible
+*
+      IF( N.EQ.0 )
+     $   RETURN
+*
+      IF( LSAME( DIRECT, 'F' ) ) THEN
+         PREVLASTV = N
+         DO I = 1, K
+            PREVLASTV = MAX( PREVLASTV, I )
+            IF( TAU( I ).EQ.ZERO ) THEN
+*
+*              H(i)  =  I
+*
+               DO J = 1, I
+                  T( J, I ) = ZERO
+               END DO
+            ELSE
+*
+*              general case
+*
+               IF( LSAME( STOREV, 'C' ) ) THEN
+*                 Skip any trailing zeros.
+                  DO LASTV = N, I+1, -1
+                     IF( V( LASTV, I ).NE.ZERO ) EXIT
+                  END DO
+                  DO J = 1, I-1
+                     T( J, I ) = -TAU( I ) * CONJG( V( I , J ) )
+                  END DO                     
+                  J = MIN( LASTV, PREVLASTV )
+*
+*                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i)
+*
+                  CALL ZGEMV( 'Conjugate transpose', J-I, I-1,
+     $                        -TAU( I ), V( I+1, 1 ), LDV, 
+     $                        V( I+1, I ), 1, ONE, T( 1, I ), 1 )
+               ELSE
+*                 Skip any trailing zeros.
+                  DO LASTV = N, I+1, -1
+                     IF( V( I, LASTV ).NE.ZERO ) EXIT
+                  END DO
+                  DO J = 1, I-1
+                     T( J, I ) = -TAU( I ) * V( J , I )
+                  END DO                     
+                  J = MIN( LASTV, PREVLASTV )
+*
+*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H
+*
+                  CALL ZGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ),
+     $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV,
+     $                        ONE, T( 1, I ), LDT )                  
+               END IF
+*
+*              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
+*
+               CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
+     $                     LDT, T( 1, I ), 1 )
+               T( I, I ) = TAU( I )
+               IF( I.GT.1 ) THEN
+                  PREVLASTV = MAX( PREVLASTV, LASTV )
+               ELSE
+                  PREVLASTV = LASTV
+               END IF
+             END IF
+         END DO
+      ELSE
+         PREVLASTV = 1
+         DO I = K, 1, -1
+            IF( TAU( I ).EQ.ZERO ) THEN
+*
+*              H(i)  =  I
+*
+               DO J = I, K
+                  T( J, I ) = ZERO
+               END DO
+            ELSE
+*
+*              general case
+*
+               IF( I.LT.K ) THEN
+                  IF( LSAME( STOREV, 'C' ) ) THEN
+*                    Skip any leading zeros.
+                     DO LASTV = 1, I-1
+                        IF( V( LASTV, I ).NE.ZERO ) EXIT
+                     END DO
+                     DO J = I+1, K
+                        T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) )
+                     END DO                        
+                     J = MAX( LASTV, PREVLASTV )
+*
+*                    T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i)
+*
+                     CALL ZGEMV( 'Conjugate transpose', N-K+I-J, K-I,
+     $                           -TAU( I ), V( J, I+1 ), LDV, V( J, I ),
+     $                           1, ONE, T( I+1, I ), 1 )
+                  ELSE
+*                    Skip any leading zeros.
+                     DO LASTV = 1, I-1
+                        IF( V( I, LASTV ).NE.ZERO ) EXIT
+                     END DO
+                     DO J = I+1, K
+                        T( J, I ) = -TAU( I ) * V( J, N-K+I )
+                     END DO                                           
+                     J = MAX( LASTV, PREVLASTV )
+*
+*                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H
+*
+                     CALL ZGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ),
+     $                           V( I+1, J ), LDV, V( I, J ), LDV,
+     $                           ONE, T( I+1, I ), LDT )                     
+                  END IF
+*
+*                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
+*
+                  CALL ZTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
+     $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
+                  IF( I.GT.1 ) THEN
+                     PREVLASTV = MIN( PREVLASTV, LASTV )
+                  ELSE
+                     PREVLASTV = LASTV
+                  END IF
+               END IF
+               T( I, I ) = TAU( I )
+            END IF
+         END DO
+      END IF
+      RETURN
+*
+*     End of ZLARFT
+*
+      END