/* -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib; on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */ #include "f2c.h" /* Table of constant values */ static integer c__1 = 1; static doublereal c_b10 = 0.; /* Subroutine */ int splicingdlarft_(char *direct, char *storev, integer *n, integer * k, doublereal *v, integer *ldv, doublereal *tau, doublereal *t, integer *ldt) { /* System generated locals */ integer t_dim1, t_offset, v_dim1, v_offset, i__1, i__2, i__3; doublereal d__1; /* Local variables */ static integer i__, j, prevlastv; static doublereal vii; extern logical splicinglsame_(char *, char *); extern /* Subroutine */ int splicingdgemv_(char *, integer *, integer *, doublereal *, doublereal *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *); static integer lastv; extern /* Subroutine */ int splicingdtrmv_(char *, char *, char *, integer *, doublereal *, integer *, doublereal *, integer *); /* -- LAPACK auxiliary routine (version 3.3.1) -- -- LAPACK is a software package provided by Univ. of Tennessee, -- -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -- April 2011 -- Purpose ======= 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 Arguments ========= DIRECT (input) 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) STOREV (input) CHARACTER*1 Specifies how the vectors which define the elementary reflectors are stored (see also Further Details): = 'C': columnwise = 'R': rowwise N (input) INTEGER The order of the block reflector H. N >= 0. K (input) INTEGER The order of the triangular factor T (= the number of elementary reflectors). K >= 1. V (input/output) DOUBLE PRECISION array, dimension (LDV,K) if STOREV = 'C' (LDV,N) if STOREV = 'R' The matrix V. See further details. LDV (input) INTEGER The leading dimension of the array V. If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. TAU (input) DOUBLE PRECISION array, dimension (K) TAU(i) must contain the scalar factor of the elementary reflector H(i). T (output) 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. LDT (input) INTEGER The leading dimension of the array T. LDT >= K. Further Details =============== 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 ) ===================================================================== Quick return if possible Parameter adjustments */ v_dim1 = *ldv; v_offset = 1 + v_dim1; v -= v_offset; --tau; t_dim1 = *ldt; t_offset = 1 + t_dim1; t -= t_offset; /* Function Body */ if (*n == 0) { return 0; } if (splicinglsame_(direct, "F")) { prevlastv = *n; i__1 = *k; for (i__ = 1; i__ <= i__1; ++i__) { prevlastv = max(i__,prevlastv); if (tau[i__] == 0.) { /* H(i) = I */ i__2 = i__; for (j = 1; j <= i__2; ++j) { t[j + i__ * t_dim1] = 0.; /* L10: */ } } else { /* general case */ vii = v[i__ + i__ * v_dim1]; v[i__ + i__ * v_dim1] = 1.; if (splicinglsame_(storev, "C")) { /* Skip any trailing zeros. */ lastv = *n; L14: if (v[lastv + i__ * v_dim1] != 0.) { goto L15; } if (lastv == i__ + 1) { goto L15; } --lastv; goto L14; L15: /* DO LASTV = N, I+1, -1 IF( V( LASTV, I ).NE.ZERO ) EXIT END DO */ j = min(lastv,prevlastv); /* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i) */ i__2 = j - i__ + 1; i__3 = i__ - 1; d__1 = -tau[i__]; splicingdgemv_("Transpose", &i__2, &i__3, &d__1, &v[i__ + v_dim1], ldv, &v[i__ + i__ * v_dim1], &c__1, &c_b10, &t[ i__ * t_dim1 + 1], &c__1); } else { /* Skip any trailing zeros. */ lastv = *n; L16: if (v[i__ + lastv * v_dim1] != 0.) { goto L17; } if (lastv == i__ + 1) { goto L17; } --lastv; goto L16; L17: /* DO LASTV = N, I+1, -1 IF( V( I, LASTV ).NE.ZERO ) EXIT END DO */ j = min(lastv,prevlastv); /* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T */ i__2 = i__ - 1; i__3 = j - i__ + 1; d__1 = -tau[i__]; splicingdgemv_("No transpose", &i__2, &i__3, &d__1, &v[i__ * v_dim1 + 1], ldv, &v[i__ + i__ * v_dim1], ldv, & c_b10, &t[i__ * t_dim1 + 1], &c__1); } v[i__ + i__ * v_dim1] = vii; /* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) */ i__2 = i__ - 1; splicingdtrmv_("Upper", "No transpose", "Non-unit", &i__2, &t[ t_offset], ldt, &t[i__ * t_dim1 + 1], &c__1); t[i__ + i__ * t_dim1] = tau[i__]; if (i__ > 1) { prevlastv = max(prevlastv,lastv); } else { prevlastv = lastv; } } /* L20: */ } } else { prevlastv = 1; for (i__ = *k; i__ >= 1; --i__) { if (tau[i__] == 0.) { /* H(i) = I */ i__1 = *k; for (j = i__; j <= i__1; ++j) { t[j + i__ * t_dim1] = 0.; /* L30: */ } } else { /* general case */ if (i__ < *k) { if (splicinglsame_(storev, "C")) { vii = v[*n - *k + i__ + i__ * v_dim1]; v[*n - *k + i__ + i__ * v_dim1] = 1.; /* Skip any leading zeros. */ lastv = 1; L34: if (v[lastv + i__ * v_dim1] != 0.) { goto L35; } if (lastv == i__ - 1) { goto L35; } ++lastv; goto L34; L35: /* DO LASTV = 1, I-1 IF( V( LASTV, I ).NE.ZERO ) EXIT 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) */ i__1 = *n - *k + i__ - j + 1; i__2 = *k - i__; d__1 = -tau[i__]; splicingdgemv_("Transpose", &i__1, &i__2, &d__1, &v[j + (i__ + 1) * v_dim1], ldv, &v[j + i__ * v_dim1], & c__1, &c_b10, &t[i__ + 1 + i__ * t_dim1], & c__1); v[*n - *k + i__ + i__ * v_dim1] = vii; } else { vii = v[i__ + (*n - *k + i__) * v_dim1]; v[i__ + (*n - *k + i__) * v_dim1] = 1.; /* Skip any leading zeros. */ lastv = 1; /* L36: */ if (v[i__ + lastv * v_dim1] != 0.) { goto L37; } if (lastv == i__ - 1) { goto L37; } ++lastv; L37: /* DO LASTV = 1, I-1 IF( V( I, LASTV ).NE.ZERO ) EXIT 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 */ i__1 = *k - i__; i__2 = *n - *k + i__ - j + 1; d__1 = -tau[i__]; splicingdgemv_("No transpose", &i__1, &i__2, &d__1, &v[i__ + 1 + j * v_dim1], ldv, &v[i__ + j * v_dim1], ldv, &c_b10, &t[i__ + 1 + i__ * t_dim1], & c__1); v[i__ + (*n - *k + i__) * v_dim1] = vii; } /* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) */ i__1 = *k - i__; splicingdtrmv_("Lower", "No transpose", "Non-unit", &i__1, &t[i__ + 1 + (i__ + 1) * t_dim1], ldt, &t[i__ + 1 + i__ * t_dim1], &c__1) ; if (i__ > 1) { prevlastv = min(prevlastv,lastv); } else { prevlastv = lastv; } } t[i__ + i__ * t_dim1] = tau[i__]; } /* L40: */ } } return 0; /* End of DLARFT */ } /* splicingdlarft_ */