1 // =============================================================================
2 // === spqr_larftb =============================================================
3 // =============================================================================
4
5 // Apply a set of Householder reflections to a matrix. Given the vectors
6 // V and coefficients Tau, construct the matrix T and then apply the updates.
7 // In MATLAB (1-based indexing), this function computes the following:
8
9 /*
10 function C = larftb (C, V, Tau, method)
11 [v k] = size (V) ;
12 [m n] = size (C) ;
13 % construct T for the compact WY representation
14 V = tril (V,-1) + eye (v,k) ;
15 T = zeros (k,k) ;
16 T (1,1) = Tau (1) ;
17 for j = 2:k
18 tau = Tau (j) ;
19 z = -tau * V (:, 1:j-1)' * V (:,j) ;
20 T (1:j-1,j) = T (1:j-1,1:j-1) * z ;
21 T (j,j) = tau ;
22 end
23 % apply the updates
24 if (method == 0)
25 C = C - V * T' * V' * C ; % method 0: Left, Transpose
26 elseif (method == 1)
27 C = C - V * T * V' * C ; % method 1: Left, No Transpose
28 elseif (method == 2)
29 C = C - C * V * T' * V' ; % method 2: Right, Transpose
30 elseif (method == 3)
31 C = C - C * V * T * V' ; % method 3: Right, No Transpose
32 end
33 */
34
35 #include "spqr.hpp"
36
spqr_private_larft(char direct,char storev,Long n,Long k,double * V,Long ldv,double * Tau,double * T,Long ldt,cholmod_common * cc)37 inline void spqr_private_larft (char direct, char storev, Long n, Long k,
38 double *V, Long ldv, double *Tau, double *T, Long ldt, cholmod_common *cc)
39 {
40 BLAS_INT N = n, K = k, LDV = ldv, LDT = ldt ;
41 if (CHECK_BLAS_INT &&
42 !(EQ (N,n) && EQ (K,k) && EQ (LDV,ldv) && EQ (LDT,ldt)))
43 {
44 cc->blas_ok = FALSE ;
45 }
46 if (!CHECK_BLAS_INT || cc->blas_ok)
47 {
48 LAPACK_DLARFT (&direct, &storev, &N, &K, V, &LDV, Tau, T, &LDT) ;
49 }
50 }
51
spqr_private_larft(char direct,char storev,Long n,Long k,Complex * V,Long ldv,Complex * Tau,Complex * T,Long ldt,cholmod_common * cc)52 inline void spqr_private_larft (char direct, char storev, Long n, Long k,
53 Complex *V, Long ldv, Complex *Tau, Complex *T, Long ldt,
54 cholmod_common *cc)
55 {
56 BLAS_INT N = n, K = k, LDV = ldv, LDT = ldt ;
57 if (CHECK_BLAS_INT &&
58 !(EQ (N,n) && EQ (K,k) && EQ (LDV,ldv) && EQ (LDT,ldt)))
59 {
60 cc->blas_ok = FALSE ;
61 }
62 if (!CHECK_BLAS_INT || cc->blas_ok)
63 {
64 LAPACK_ZLARFT (&direct, &storev, &N, &K, V, &LDV, Tau, T, &LDT) ;
65 }
66 }
67
68
spqr_private_larfb(char side,char trans,char direct,char storev,Long m,Long n,Long k,double * V,Long ldv,double * T,Long ldt,double * C,Long ldc,double * Work,Long ldwork,cholmod_common * cc)69 inline void spqr_private_larfb (char side, char trans, char direct, char storev,
70 Long m, Long n, Long k, double *V, Long ldv, double *T, Long ldt, double *C,
71 Long ldc, double *Work, Long ldwork, cholmod_common *cc)
72 {
73 BLAS_INT M = m, N = n, K = k, LDV = ldv, LDT = ldt, LDC = ldc,
74 LDWORK = ldwork ;
75 if (CHECK_BLAS_INT &&
76 !(EQ (M,m) && EQ (N,n) && EQ (K,k) && EQ (LDV,ldv) &&
77 EQ (LDT,ldt) && EQ (LDV,ldv) && EQ (LDWORK,ldwork)))
78 {
79 cc->blas_ok = FALSE ;
80 }
81 if (!CHECK_BLAS_INT || cc->blas_ok)
82 {
83 LAPACK_DLARFB (&side, &trans, &direct, &storev, &M, &N, &K, V, &LDV,
84 T, &LDT, C, &LDC, Work, &LDWORK) ;
85 }
86 }
87
88
spqr_private_larfb(char side,char trans,char direct,char storev,Long m,Long n,Long k,Complex * V,Long ldv,Complex * T,Long ldt,Complex * C,Long ldc,Complex * Work,Long ldwork,cholmod_common * cc)89 inline void spqr_private_larfb (char side, char trans, char direct, char storev,
90 Long m, Long n, Long k, Complex *V, Long ldv, Complex *T, Long ldt,
91 Complex *C, Long ldc, Complex *Work, Long ldwork, cholmod_common *cc)
92 {
93 char tr = (trans == 'T') ? 'C' : 'N' ; // change T to C
94 BLAS_INT M = m, N = n, K = k, LDV = ldv, LDT = ldt, LDC = ldc,
95 LDWORK = ldwork ;
96 if (CHECK_BLAS_INT &&
97 !(EQ (M,m) && EQ (N,n) && EQ (K,k) && EQ (LDV,ldv) &&
98 EQ (LDT,ldt) && EQ (LDV,ldv) && EQ (LDWORK,ldwork)))
99 {
100 cc->blas_ok = FALSE ;
101 }
102 if (!CHECK_BLAS_INT || cc->blas_ok)
103 {
104 LAPACK_ZLARFB (&side, &tr, &direct, &storev, &M, &N, &K, V, &LDV,
105 T, &LDT, C, &LDC, Work, &LDWORK) ;
106 }
107 }
108
109
110 // =============================================================================
111
spqr_larftb(int method,Long m,Long n,Long k,Long ldc,Long ldv,Entry * V,Entry * Tau,Entry * C,Entry * W,cholmod_common * cc)112 template <typename Entry> void spqr_larftb
113 (
114 // inputs, not modified (V is modified and then restored on output)
115 int method, // 0,1,2,3
116 Long m, // C is m-by-n
117 Long n,
118 Long k, // V is v-by-k
119 // for methods 0 and 1, v = m,
120 // for methods 2 and 3, v = n
121 Long ldc, // leading dimension of C
122 Long ldv, // leading dimension of V
123 Entry *V, // V is v-by-k, unit lower triangular (diag not stored)
124 Entry *Tau, // size k, the k Householder coefficients
125
126 // input/output
127 Entry *C, // C is m-by-n, with leading dimension ldc
128
129 // workspace, not defined on input or output
130 Entry *W, // for methods 0,1: size k*k + n*k
131 // for methods 2,3: size k*k + m*k
132 cholmod_common *cc
133 )
134 {
135 Entry *T, *Work ;
136
137 // -------------------------------------------------------------------------
138 // check inputs and split up workspace
139 // -------------------------------------------------------------------------
140
141 if (m <= 0 || n <= 0 || k <= 0)
142 {
143 return ; // nothing to do
144 }
145
146 T = W ; // triangular k-by-k matrix for block reflector
147 Work = W + k*k ; // workspace of size n*k or m*k for larfb
148
149 // -------------------------------------------------------------------------
150 // construct and apply the k-by-k upper triangular matrix T
151 // -------------------------------------------------------------------------
152
153 // larft and larfb are always used "Forward" and "Columnwise"
154
155 if (method == SPQR_QTX)
156 {
157 ASSERT (m >= k) ;
158 spqr_private_larft ('F', 'C', m, k, V, ldv, Tau, T, k, cc) ;
159 // Left, Transpose, Forward, Columwise:
160 spqr_private_larfb ('L', 'T', 'F', 'C', m, n, k, V, ldv, T, k, C, ldc,
161 Work, n, cc) ;
162 }
163 else if (method == SPQR_QX)
164 {
165 ASSERT (m >= k) ;
166 spqr_private_larft ('F', 'C', m, k, V, ldv, Tau, T, k, cc) ;
167 // Left, No Transpose, Forward, Columwise:
168 spqr_private_larfb ('L', 'N', 'F', 'C', m, n, k, V, ldv, T, k, C, ldc,
169 Work, n, cc) ;
170 }
171 else if (method == SPQR_XQT)
172 {
173 ASSERT (n >= k) ;
174 spqr_private_larft ('F', 'C', n, k, V, ldv, Tau, T, k, cc) ;
175 // Right, Transpose, Forward, Columwise:
176 spqr_private_larfb ('R', 'T', 'F', 'C', m, n, k, V, ldv, T, k, C, ldc,
177 Work, m, cc) ;
178 }
179 else if (method == SPQR_XQ)
180 {
181 ASSERT (n >= k) ;
182 spqr_private_larft ('F', 'C', n, k, V, ldv, Tau, T, k, cc) ;
183 // Right, No Transpose, Forward, Columwise:
184 spqr_private_larfb ('R', 'N', 'F', 'C', m, n, k, V, ldv, T, k, C, ldc,
185 Work, m, cc) ;
186 }
187 }
188
189 // =============================================================================
190
191 template void spqr_larftb <double>
192 (
193 // inputs, not modified (V is modified and then restored on output)
194 int method, // 0,1,2,3
195 Long m, // C is m-by-n
196 Long n,
197 Long k, // V is v-by-k
198 // for methods 0 and 1, v = m,
199 // for methods 2 and 3, v = n
200 Long ldc, // leading dimension of C
201 Long ldv, // leading dimension of V
202 double *V, // V is v-by-k, unit lower triangular (diag not stored)
203 double *Tau, // size k, the k Householder coefficients
204
205 // input/output
206 double *C, // C is m-by-n, with leading dimension ldc
207
208 // workspace, not defined on input or output
209 double *W, // for methods 0,1: size k*k + n*k
210 // for methods 2,3: size k*k + m*k
211 cholmod_common *cc
212 ) ;
213
214 // =============================================================================
215
216 template void spqr_larftb <Complex>
217 (
218 // inputs, not modified (V is modified and then restored on output)
219 int method, // 0,1,2,3
220 Long m, // C is m-by-n
221 Long n,
222 Long k, // V is v-by-k
223 // for methods 0 and 1, v = m,
224 // for methods 2 and 3, v = n
225 Long ldc, // leading dimension of C
226 Long ldv, // leading dimension of V
227 Complex *V, // V is v-by-k, unit lower triangular (diag not stored)
228 Complex *Tau, // size k, the k Householder coefficients
229
230 // input/output
231 Complex *C, // C is m-by-n, with leading dimension ldc
232
233 // workspace, not defined on input or output
234 Complex *W, // for methods 0,1: size k*k + n*k
235 // for methods 2,3: size k*k + m*k
236 cholmod_common *cc
237 ) ;
238