Line data Source code
1 : /*
2 : * The MIT License (MIT)
3 : *
4 : * Copyright (c) 2017 Pantelis Sopasakis (https://alphaville.github.io),
5 : * Krina Menounou (https://www.linkedin.com/in/krinamenounou),
6 : * Panagiotis Patrinos (http://homes.esat.kuleuven.be/~ppatrino)
7 : *
8 : * Permission is hereby granted, free of charge, to any person obtaining a copy
9 : * of this software and associated documentation files (the "Software"), to deal
10 : * in the Software without restriction, including without limitation the rights
11 : * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 : * copies of the Software, and to permit persons to whom the Software is
13 : * furnished to do so, subject to the following conditions:
14 : *
15 : * The above copyright notice and this permission notice shall be included in all
16 : * copies or substantial portions of the Software.
17 : *
18 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 : * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 : * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 : * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24 : * SOFTWARE.
25 : *
26 : */
27 : #include "directions.h"
28 :
29 :
30 : /*TODO Allocate this variable within work */
31 : static scs_float * HY; /* Vector H*Y of length l */
32 :
33 3620 : scs_int scs_reset_direction_cache(ScsDirectionCache * cache) {
34 3620 : cache->mem_cursor = 0; /* set active memory to 0 */
35 3620 : cache->current_mem = 0;
36 3620 : return SCS_DIRECTION_CACHE_RESET;
37 : }
38 :
39 41 : scs_int scs_compute_dir_anderson(ScsWork *work) {
40 : /* --- DECLARATIONS --- */
41 : ScsDirectionCache * cache; /* the SU cache (pointer) */
42 41 : const scs_int l = work->l; /* size of vectors */
43 : scs_float * s_current;
44 : scs_float * y_current;
45 : scs_float * s_minus_y_current;
46 : scs_int colsY;
47 :
48 : #ifdef SVD_ACTIVATED
49 41 : scs_float rcond = 1e-8;
50 : scs_int rank;
51 : scs_float * singular_values;
52 : #endif
53 :
54 : scs_float * copy_cache_Y;
55 :
56 41 : cache = work->direction_cache;
57 :
58 : #ifdef SVD_ACTIVATED
59 41 : singular_values = cache->ls_wspace + cache->ls_wspace_length;
60 41 : copy_cache_Y = singular_values + cache->mem;
61 : #else
62 : copy_cache_Y = cache->ls_wspace + cache->ls_wspace_length;
63 : #endif
64 :
65 :
66 :
67 : /* d [work->dir] = -R [work->R] */
68 41 : scs_set_as_scaled_array(work->dir, work->R, -1.0, l);
69 :
70 41 : s_current = cache->S + (cache->mem_cursor * l);
71 41 : y_current = cache->U + (cache->mem_cursor * l);
72 41 : s_minus_y_current = cache->S_minus_Y + (cache->mem_cursor * l);
73 :
74 41 : memcpy(s_current, work->Sk, l * sizeof (scs_float));
75 41 : memcpy(y_current, work->Yk, l * sizeof (scs_float));
76 :
77 41 : cache->current_mem++; /* increment the current memory */
78 41 : if (cache->current_mem > cache->mem)
79 18 : cache->current_mem = cache->mem;
80 :
81 41 : colsY = cache->current_mem;
82 :
83 : /* Construct matrix S-Y by subtracting one pair of vecs at a time */
84 41 : scs_axpy(s_minus_y_current, s_current, y_current, 1.0, -1.0, l);
85 :
86 : /* Solve Yt = R with SVD (we need to copy Y into a different memory positions
87 : * because svdls modifies it (see the documentation of svdls). */
88 41 : memcpy(cache->t, work->R, l * sizeof (scs_float));
89 41 : memcpy(copy_cache_Y, cache->U, l * colsY * sizeof (scs_float));
90 :
91 : #ifdef SVD_ACTIVATED
92 : /* Solve LS with SVD */
93 41 : scs_svdls(l, colsY,
94 : copy_cache_Y,
95 : cache->t,
96 : cache->ls_wspace,
97 : cache->ls_wspace_length,
98 : rcond,
99 : singular_values,
100 : &rank);
101 : #else
102 : /* Solve LS with QR */
103 : scs_qrls(l, colsY,
104 : copy_cache_Y,
105 : cache->t,
106 : cache->ls_wspace,
107 : cache->ls_wspace_length);
108 : #endif
109 :
110 : /* dir = dir - S_minus_Y * t */
111 82 : scs_matrix_multiply(l, 1, colsY, -1.0,
112 41 : cache->S_minus_Y, 1.0, cache->t, work->dir);
113 :
114 41 : cache->mem_cursor++; /* move the cursor */
115 41 : if (cache->mem_cursor >= cache->mem)
116 7 : cache->mem_cursor = 0;
117 :
118 41 : return SCS_DIRECTION_SUCCESS;
119 : }
120 :
121 34499 : scs_int scs_compute_dir_restarted_broyden(ScsWork *work) {
122 : /* --- DECLARATIONS --- */
123 : ScsDirectionCache * cache; /* the SU cache (pointer) */
124 : scs_int i; /* index */
125 : scs_float * s_tilde_current; /* s_tilde (which is updated) */
126 : scs_float * u_new; /* new value of u */
127 : scs_float ip; /* temporary float to store inner products */
128 : scs_float s_norm_sq; /* scalar gamma as in (6.5e) */
129 34499 : scs_float theta = 0; /* theta */
130 34499 : const scs_int l = work->l; /* size of vectors */
131 34499 : const scs_float theta_bar = work->stgs->thetabar; /* parameter in Powell's trick */
132 :
133 34499 : cache = work->direction_cache; /* cache of Sk and Uk */
134 :
135 : /* d [work->dir] = -R [work->R] */
136 34499 : scs_set_as_scaled_array(work->dir, work->R, -1.0, l);
137 :
138 : /* s_tilde_current = y [work->Yk] */
139 : /* use the end of the cache to store s_tilde_current */
140 : /* later we use the same position of the S-buffer to store the current Sk */
141 34499 : s_tilde_current = cache->S + (cache->mem_cursor * l);
142 34499 : memcpy(s_tilde_current, work->Yk, l * sizeof (scs_float));
143 :
144 : /* update s and d */
145 188240 : for (i = 0; i < cache->mem_cursor; ++i) {
146 : scs_float * s_i; /* pointer to the current value of s_i */
147 : scs_float * u_i; /* pointer to the current value of u_i */
148 153741 : s_i = cache->S + i * l; /* retrieve s_i from the cache */
149 153741 : u_i = cache->U + i * l; /* retrieve u_i from the cache */
150 153741 : ip = scs_inner_product(s_i, s_tilde_current, l);
151 153741 : scs_add_scaled_array(s_tilde_current, u_i, l, ip); /* update s_tilde */
152 153741 : ip = scs_inner_product(s_i, work->dir, l);
153 153741 : scs_add_scaled_array(work->dir, u_i, l, ip); /* update direction */
154 : }
155 :
156 : /* compute theta */
157 34499 : ip = scs_inner_product(s_tilde_current, work->Sk, l);
158 34499 : s_norm_sq = scs_norm_squared(work->Sk, l);
159 :
160 34499 : if (ABS(ip) >= theta_bar * s_norm_sq) {
161 : theta = 1;
162 : } else {
163 3113 : theta = s_norm_sq * (1 - SGN(ip) * theta_bar) / (s_norm_sq - ip);
164 : /* s_tilde_current = (1-theta)*s + theta*s_tilde_current */
165 3113 : scs_axpy(s_tilde_current, s_tilde_current, work->Sk, theta, 1 - theta, l);
166 : }
167 :
168 : /* FINALISE */
169 :
170 : /* update u_new (at the end of the buffer) */
171 34499 : u_new = cache->U + (cache->mem_cursor * l);
172 34499 : ip = (1 - theta) * s_norm_sq + theta * ip;
173 377589 : for (i = 0; i < l; ++i) {
174 343090 : u_new[i] = (work->Sk[i] - s_tilde_current[i]) / ip;
175 : }
176 : /* update direction */
177 34499 : ip = scs_inner_product(work->Sk, work->dir, l); /* s'd */
178 34499 : scs_add_scaled_array(work->dir, u_new, l, ip);
179 :
180 : /* push s into the buffer */
181 34499 : memcpy(s_tilde_current, work->Sk, l * sizeof (scs_float));
182 :
183 34499 : cache->mem_cursor++; /* move the cursor */
184 :
185 : /* if the cursor has exceeded the last position, reset the cache */
186 34499 : if (cache->mem_cursor >= cache->mem) {
187 3295 : return scs_reset_direction_cache(cache); /* returns SCS_DIRECTION_CACHE_RESET */
188 : }
189 :
190 : return SCS_DIRECTION_CACHE_INCREMENT;
191 : }
192 :
193 : /* LCOV_EXCL_START */
194 : scs_int scs_compute_dir_full_broyden(ScsWork *work, scs_int i) {
195 : scs_float ip = 0;
196 : scs_float tmp = 0;
197 :
198 : if (i == 0 || HY == SCS_NULL) {
199 : /* HY is allocated the first time this function is called (that is, for i==0) */
200 : HY = malloc(work->l * sizeof (*HY));
201 : if (HY == SCS_NULL) {
202 : scs_printf("ERROR: allocating `HY` in `scs_compute_dir_full_broyden` failure\n");
203 : return SCS_DIRECTION_ERROR;
204 : }
205 : }
206 :
207 : if ((work->stgs->broyden_init_scaling && i == 1)
208 : || (work->stgs->tRule == 1 || work->stgs->tRule == 2)) {
209 : ip = scs_inner_product(work->Yk, work->Sk, work->l);
210 : }
211 :
212 : if (work->stgs->broyden_init_scaling && i == 1) {
213 : scs_int i;
214 : tmp = ip / scs_norm(work->Yk, work->l);
215 : for (i = 0; i < work->l; ++i) {
216 : work->H[i * (work->l + 1)] = tmp;
217 : }
218 : }
219 :
220 : return SCS_DIRECTION_SUCCESS;
221 : }
222 :
223 : /* LCOV_EXCL_STOP */
224 :
225 29120 : scs_int scs_compute_direction(ScsWork *work, scs_int i) {
226 29120 : scs_int status = SCS_DIRECTION_SUCCESS;
227 :
228 29120 : switch (work->stgs->direction) {
229 : case fixed_point_residual:
230 89 : scs_set_as_scaled_array(work->dir, work->R, -1.0, work->l); /* dir = -R */
231 89 : status = SCS_DIRECTION_SUCCESS;
232 89 : break;
233 : case restarted_broyden:
234 28990 : status = scs_compute_dir_restarted_broyden(work);
235 28990 : break;
236 : case anderson_acceleration:
237 41 : status = scs_compute_dir_anderson(work);
238 41 : break;
239 : case full_broyden:
240 0 : status = scs_compute_dir_full_broyden(work, i);
241 0 : break;
242 : default:
243 : /* Not implemented yet */
244 : status = SCS_DIRECTION_ERROR;
245 : }
246 :
247 29120 : return status;
248 : }
249 :
250 : /* LCOV_EXCL_START */
251 : void scs_free_full_broyden() {
252 : if (HY != SCS_NULL) {
253 : scs_free(HY);
254 : }
255 : }
256 : /* LCOV_EXCL_STOP */
|
