1 #include "cs.h"
2 /* create a Householder reflection [v,beta,s]=house(x), overwrite x with v,
3 * where (I-beta*v*v')*x = s*e1 and e1 = [1 0 ... 0]'.
4 * Note that this CXSparse version is different than CSparse. See Higham,
5 * Accuracy & Stability of Num Algorithms, 2nd ed, 2002, page 357. */
cs_house(CS_ENTRY * x,double * beta,CS_INT n)6 CS_ENTRY cs_house (CS_ENTRY *x, double *beta, CS_INT n)
7 {
8 CS_ENTRY s = 0 ;
9 CS_INT i ;
10 if (!x || !beta) return (-1) ; /* check inputs */
11 /* s = norm(x) */
12 for (i = 0 ; i < n ; i++) s += x [i] * CS_CONJ (x [i]) ;
13 s = sqrt (s) ;
14 if (s == 0)
15 {
16 (*beta) = 0 ;
17 x [0] = 1 ;
18 }
19 else
20 {
21 /* s = sign(x[0]) * norm (x) ; */
22 if (x [0] != 0)
23 {
24 s *= x [0] / CS_ABS (x [0]) ;
25 }
26 x [0] += s ;
27 (*beta) = 1. / CS_REAL (CS_CONJ (s) * x [0]) ;
28 }
29 return (-s) ;
30 }
31