1 // =============================================================================
2 // === spqr_happly_work ========================================================
3 // =============================================================================
4
5 // Determines the workspace workspace needed by spqr-happly
6
7 #include "spqr.hpp"
8
spqr_happly_work(int method,Long m,Long n,Long nh,Long * Hp,Long hchunk,Long * p_vmax,Long * p_vsize,Long * p_csize)9 int spqr_happly_work
10 (
11 // input
12 int method, // 0,1,2,3
13
14 Long m, // X is m-by-n
15 Long n,
16
17 // FUTURE : make H cholmod_sparse:
18 Long nh, // number of Householder vectors
19 Long *Hp, // size nh+1, column pointers for H
20 Long hchunk,
21
22 // outputs; sizes of workspaces needed
23 Long *p_vmax,
24 Long *p_vsize,
25 Long *p_csize
26 )
27 {
28 Long maxhlen, h, hlen, vmax, mh, vsize, csize, vsize1, vsize2 ;
29 int ok = TRUE ;
30
31 // -------------------------------------------------------------------------
32 // get inputs
33 // -------------------------------------------------------------------------
34
35 *p_vmax = 0 ;
36 *p_vsize = 0 ;
37 *p_csize = 0 ;
38
39 if (m == 0 || n == 0 || nh == 0)
40 {
41 // nothing to do
42 return (TRUE) ;
43 }
44
45 // -------------------------------------------------------------------------
46 // determine the length of the longest Householder vector
47 // -------------------------------------------------------------------------
48
49 maxhlen = 1 ;
50 for (h = 0 ; h < nh ; h++)
51 {
52 hlen = Hp [h+1] - Hp [h] ;
53 maxhlen = MAX (maxhlen, hlen) ;
54 }
55
56 // number of rows of H
57 mh = (method == 0 || method == 1) ? m : n ;
58
59 // -------------------------------------------------------------------------
60 // determine workspace sizes
61 // -------------------------------------------------------------------------
62
63 // Long overflow cannot occur with vmax since H is already allocated
64 if (method == 0 || method == 3)
65 {
66 // apply H in the forward direction; H(0) first, H(nh-1) last
67 vmax = 2 * maxhlen + 8 ;
68 }
69 else
70 {
71 // apply H in the backward direction; H(nh-1) first, H(0) last
72 vmax = maxhlen + hchunk ;
73 }
74
75 vmax = MIN (vmax, mh) ;
76 vmax = MAX (vmax, 2) ;
77
78 // csize = vmax * ((method <= 1) ? n : m) ;
79 csize = spqr_mult (vmax, (method <= 1) ? n : m, &ok) ;
80
81 // vsize = (hchunk*hchunk + ((method <= 1) ? n : m)*hchunk + vmax*hchunk) ;
82 vsize = spqr_mult (hchunk, hchunk, &ok) ;
83 vsize1 = spqr_mult ((method <= 1) ? n : m, hchunk, &ok) ;
84 vsize2 = spqr_mult (vmax, hchunk, &ok) ;
85 vsize = spqr_add (vsize, vsize1, &ok) ;
86 vsize = spqr_add (vsize, vsize2, &ok) ;
87
88 // -------------------------------------------------------------------------
89 // return workspace sizes
90 // -------------------------------------------------------------------------
91
92 *p_vmax = vmax ;
93 *p_vsize = vsize ;
94 *p_csize = csize ;
95 return (ok) ;
96 }
97