X-Git-Url: https://vcs.maemo.org/git/?a=blobdiff_plain;f=3rdparty%2Flapack%2Fdlarft.c;fp=3rdparty%2Flapack%2Fdlarft.c;h=62f8753a2601bb55a96e412498f7c57df2a63d85;hb=e4c14cdbdf2fe805e79cd96ded236f57e7b89060;hp=0000000000000000000000000000000000000000;hpb=454138ff8a20f6edb9b65a910101403d8b520643;p=opencv diff --git a/3rdparty/lapack/dlarft.c b/3rdparty/lapack/dlarft.c new file mode 100644 index 0000000..62f8753 --- /dev/null +++ b/3rdparty/lapack/dlarft.c @@ -0,0 +1,266 @@ +#include "clapack.h" + +/* Table of constant values */ + +static integer c__1 = 1; +static doublereal c_b8 = 0.; + +/* Subroutine */ int dlarft_(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 */ + integer i__, j; + doublereal vii; + extern logical lsame_(char *, char *); + extern /* Subroutine */ int dgemv_(char *, integer *, integer *, + doublereal *, doublereal *, integer *, doublereal *, integer *, + doublereal *, doublereal *, integer *), dtrmv_(char *, + char *, char *, integer *, doublereal *, integer *, doublereal *, + integer *); + + +/* -- LAPACK auxiliary routine (version 3.1) -- */ +/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ +/* November 2006 */ + +/* .. Scalar Arguments .. */ +/* .. */ +/* .. Array Arguments .. */ +/* .. */ + +/* 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' */ + +/* 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 * 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 ) */ + +/* ===================================================================== */ + +/* .. Parameters .. */ +/* .. */ +/* .. Local Scalars .. */ +/* .. */ +/* .. External Subroutines .. */ +/* .. */ +/* .. External Functions .. */ +/* .. */ +/* .. Executable Statements .. */ + +/* 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 (lsame_(direct, "F")) { + i__1 = *k; + for (i__ = 1; i__ <= i__1; ++i__) { + 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 (lsame_(storev, "C")) { + +/* T(1:i-1,i) := - tau(i) * V(i:n,1:i-1)' * V(i:n,i) */ + + i__2 = *n - i__ + 1; + i__3 = i__ - 1; + d__1 = -tau[i__]; + dgemv_("Transpose", &i__2, &i__3, &d__1, &v[i__ + v_dim1], + ldv, &v[i__ + i__ * v_dim1], &c__1, &c_b8, &t[ + i__ * t_dim1 + 1], &c__1); + } else { + +/* T(1:i-1,i) := - tau(i) * V(1:i-1,i:n) * V(i,i:n)' */ + + i__2 = i__ - 1; + i__3 = *n - i__ + 1; + d__1 = -tau[i__]; + dgemv_("No transpose", &i__2, &i__3, &d__1, &v[i__ * + v_dim1 + 1], ldv, &v[i__ + i__ * v_dim1], ldv, & + c_b8, &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; + dtrmv_("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__]; + } +/* L20: */ + } + } else { + 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 (lsame_(storev, "C")) { + vii = v[*n - *k + i__ + i__ * v_dim1]; + v[*n - *k + i__ + i__ * v_dim1] = 1.; + +/* T(i+1:k,i) := */ +/* - tau(i) * V(1:n-k+i,i+1:k)' * V(1:n-k+i,i) */ + + i__1 = *n - *k + i__; + i__2 = *k - i__; + d__1 = -tau[i__]; + dgemv_("Transpose", &i__1, &i__2, &d__1, &v[(i__ + 1) + * v_dim1 + 1], ldv, &v[i__ * v_dim1 + 1], & + c__1, &c_b8, &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.; + +/* T(i+1:k,i) := */ +/* - tau(i) * V(i+1:k,1:n-k+i) * V(i,1:n-k+i)' */ + + i__1 = *k - i__; + i__2 = *n - *k + i__; + d__1 = -tau[i__]; + dgemv_("No transpose", &i__1, &i__2, &d__1, &v[i__ + + 1 + v_dim1], ldv, &v[i__ + v_dim1], ldv, & + c_b8, &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__; + dtrmv_("Lower", "No transpose", "Non-unit", &i__1, &t[i__ + + 1 + (i__ + 1) * t_dim1], ldt, &t[i__ + 1 + i__ * + t_dim1], &c__1) + ; + } + t[i__ + i__ * t_dim1] = tau[i__]; + } +/* L40: */ + } + } + return 0; + +/* End of DLARFT */ + +} /* dlarft_ */