SHOGUN
v1.1.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Purpose: solves quadratic programming problem for pattern recognition 00008 * for support vectors 00009 * 00010 * Written (W) 1997-1998 Alex J. Smola 00011 * Written (W) 1999-2009 Soeren Sonnenburg 00012 * Written (W) 1999-2008 Gunnar Raetsch 00013 * Copyright (C) 1997-2009 Fraunhofer Institute FIRST and Max-Planck-Society 00014 */ 00015 00016 #include <shogun/lib/common.h> 00017 #include <shogun/io/SGIO.h> 00018 #include <shogun/mathematics/lapack.h> 00019 #include <shogun/mathematics/Math.h> 00020 #include <shogun/classifier/svm/pr_loqo.h> 00021 00022 #include <math.h> 00023 #include <time.h> 00024 #include <stdlib.h> 00025 #include <stdio.h> 00026 #ifdef HAVE_LAPACK 00027 extern "C" { 00028 #ifdef HAVE_ATLAS 00029 #include <atlas/clapack.h> 00030 #endif 00031 } 00032 #endif 00033 00034 namespace shogun 00035 { 00036 00037 #define PREDICTOR 1 00038 #define CORRECTOR 2 00039 00040 /***************************************************************** 00041 replace this by any other function that will exit gracefully 00042 in a larger system 00043 ***************************************************************/ 00044 00045 void nrerror(char error_text[]) 00046 { 00047 SG_SDEBUG("terminating optimizer - %s\n", error_text); 00048 // exit(1); 00049 } 00050 00051 /***************************************************************** 00052 taken from numerical recipes and modified to accept pointers 00053 moreover numerical recipes code seems to be buggy (at least the 00054 ones on the web) 00055 00056 cholesky solver and backsubstitution 00057 leaves upper right triangle intact (rows first order) 00058 ***************************************************************/ 00059 00060 #ifdef HAVE_LAPACK 00061 bool choldc(float64_t* a, int32_t n, float64_t* p) 00062 { 00063 if (n<=0) 00064 return false; 00065 00066 float64_t* a2=SG_MALLOC(float64_t, n*n); 00067 00068 for (int32_t i=0; i<n; i++) 00069 { 00070 for (int32_t j=0; j<n; j++) 00071 a2[n*i+j]=a[n*i+j]; 00072 } 00073 00074 /* int for calling external lib */ 00075 int result=clapack_dpotrf(CblasRowMajor, CblasUpper, (int) n, a2, (int) n); 00076 00077 for (int32_t i=0; i<n; i++) 00078 p[i]=a2[(n+1)*i]; 00079 00080 for (int32_t i=0; i<n; i++) 00081 { 00082 for (int32_t j=i+1; j<n; j++) 00083 { 00084 a[n*j+i]=a2[n*i+j]; 00085 } 00086 } 00087 00088 if (result>0) 00089 SG_SDEBUG("Choldc failed, matrix not positive definite\n"); 00090 00091 SG_FREE(a2); 00092 00093 return result==0; 00094 } 00095 #else 00096 bool choldc(float64_t a[], int32_t n, float64_t p[]) 00097 { 00098 void nrerror(char error_text[]); 00099 int32_t i, j, k; 00100 float64_t sum; 00101 00102 for (i = 0; i < n; i++) 00103 { 00104 for (j = i; j < n; j++) 00105 { 00106 sum=a[n*i + j]; 00107 00108 for (k=i-1; k>=0; k--) 00109 sum -= a[n*i + k]*a[n*j + k]; 00110 00111 if (i == j) 00112 { 00113 if (sum <= 0.0) 00114 { 00115 SG_SDEBUG("Choldc failed, matrix not positive definite"); 00116 sum = 0.0; 00117 return false; 00118 } 00119 00120 p[i]=sqrt(sum); 00121 00122 } 00123 else 00124 a[n*j + i] = sum/p[i]; 00125 } 00126 } 00127 00128 return true; 00129 } 00130 #endif 00131 00132 void cholsb( 00133 float64_t a[], int32_t n, float64_t p[], float64_t b[], float64_t x[]) 00134 { 00135 int32_t i, k; 00136 float64_t sum; 00137 00138 for (i=0; i<n; i++) { 00139 sum=b[i]; 00140 for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k]; 00141 x[i]=sum/p[i]; 00142 } 00143 00144 for (i=n-1; i>=0; i--) { 00145 sum=x[i]; 00146 for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k]; 00147 x[i]=sum/p[i]; 00148 } 00149 } 00150 00151 /***************************************************************** 00152 sometimes we only need the forward or backward pass of the 00153 backsubstitution, hence we provide these two routines separately 00154 ***************************************************************/ 00155 00156 void chol_forward( 00157 float64_t a[], int32_t n, float64_t p[], float64_t b[], float64_t x[]) 00158 { 00159 int32_t i, k; 00160 float64_t sum; 00161 00162 for (i=0; i<n; i++) { 00163 sum=b[i]; 00164 for (k=i-1; k>=0; k--) sum -= a[n*i + k]*x[k]; 00165 x[i]=sum/p[i]; 00166 } 00167 } 00168 00169 void chol_backward( 00170 float64_t a[], int32_t n, float64_t p[], float64_t b[], float64_t x[]) 00171 { 00172 int32_t i, k; 00173 float64_t sum; 00174 00175 for (i=n-1; i>=0; i--) { 00176 sum=b[i]; 00177 for (k=i+1; k<n; k++) sum -= a[n*k + i]*x[k]; 00178 x[i]=sum/p[i]; 00179 } 00180 } 00181 00182 /***************************************************************** 00183 solves the system | -H_x A' | |x_x| = |c_x| 00184 | A H_y| |x_y| |c_y| 00185 00186 with H_x (and H_y) positive (semidefinite) matrices 00187 and n, m the respective sizes of H_x and H_y 00188 00189 for variables see pg. 48 of notebook or do the calculations on a 00190 sheet of paper again 00191 00192 predictor solves the whole thing, corrector assues that H_x didn't 00193 change and relies on the results of the predictor. therefore do 00194 _not_ modify workspace 00195 00196 if you want to speed tune anything in the code here's the right 00197 place to do so: about 95% of the time is being spent in 00198 here. something like an iterative refinement would be nice, 00199 especially when switching from float64_t to single precision. if you 00200 have a fast parallel cholesky use it instead of the numrec 00201 implementations. 00202 00203 side effects: changes H_y (but this is just the unit matrix or zero anyway 00204 in our case) 00205 ***************************************************************/ 00206 00207 bool solve_reduced( 00208 int32_t n, int32_t m, float64_t h_x[], float64_t h_y[], float64_t a[], 00209 float64_t x_x[], float64_t x_y[], float64_t c_x[], float64_t c_y[], 00210 float64_t workspace[], int32_t step) 00211 { 00212 int32_t i,j,k; 00213 00214 float64_t *p_x; 00215 float64_t *p_y; 00216 float64_t *t_a; 00217 float64_t *t_c; 00218 float64_t *t_y; 00219 00220 p_x = workspace; /* together n + m + n*m + n + m = n*(m+2)+2*m */ 00221 p_y = p_x + n; 00222 t_a = p_y + m; 00223 t_c = t_a + n*m; 00224 t_y = t_c + n; 00225 00226 if (step == PREDICTOR) { 00227 if (!choldc(h_x, n, p_x)) /* do cholesky decomposition */ 00228 return false; 00229 00230 for (i=0; i<m; i++) /* forward pass for A' */ 00231 chol_forward(h_x, n, p_x, a+i*n, t_a+i*n); 00232 00233 for (i=0; i<m; i++) /* compute (h_y + a h_x^-1A') */ 00234 for (j=i; j<m; j++) 00235 for (k=0; k<n; k++) 00236 h_y[m*i + j] += t_a[n*j + k] * t_a[n*i + k]; 00237 00238 choldc(h_y, m, p_y); /* and cholesky decomposition */ 00239 } 00240 00241 chol_forward(h_x, n, p_x, c_x, t_c); 00242 /* forward pass for c */ 00243 00244 for (i=0; i<m; i++) { /* and solve for x_y */ 00245 t_y[i] = c_y[i]; 00246 for (j=0; j<n; j++) 00247 t_y[i] += t_a[i*n + j] * t_c[j]; 00248 } 00249 00250 cholsb(h_y, m, p_y, t_y, x_y); 00251 00252 for (i=0; i<n; i++) { /* finally solve for x_x */ 00253 t_c[i] = -t_c[i]; 00254 for (j=0; j<m; j++) 00255 t_c[i] += t_a[j*n + i] * x_y[j]; 00256 } 00257 00258 chol_backward(h_x, n, p_x, t_c, x_x); 00259 return true; 00260 } 00261 00262 /***************************************************************** 00263 matrix vector multiplication (symmetric matrix but only one triangle 00264 given). computes m*x = y 00265 no need to tune it as it's only of O(n^2) but cholesky is of 00266 O(n^3). so don't waste your time _here_ although it isn't very 00267 elegant. 00268 ***************************************************************/ 00269 00270 void matrix_vector(int32_t n, float64_t m[], float64_t x[], float64_t y[]) 00271 { 00272 int32_t i, j; 00273 00274 for (i=0; i<n; i++) { 00275 y[i] = m[(n+1) * i] * x[i]; 00276 00277 for (j=0; j<i; j++) 00278 y[i] += m[i + n*j] * x[j]; 00279 00280 for (j=i+1; j<n; j++) 00281 y[i] += m[n*i + j] * x[j]; 00282 } 00283 } 00284 00285 /***************************************************************** 00286 call only this routine; this is the only one you're interested in 00287 for doing quadratical optimization 00288 00289 the restart feature exists but it may not be of much use due to the 00290 fact that an initial setting, although close but not very close the 00291 the actual solution will result in very good starting diagnostics 00292 (primal and dual feasibility and small infeasibility gap) but incur 00293 later stalling of the optimizer afterwards as we have to enforce 00294 positivity of the slacks. 00295 ***************************************************************/ 00296 00297 int32_t pr_loqo( 00298 int32_t n, int32_t m, float64_t c[], float64_t h_x[], float64_t a[], 00299 float64_t b[], float64_t l[], float64_t u[], float64_t primal[], 00300 float64_t dual[], int32_t verb, float64_t sigfig_max, int32_t counter_max, 00301 float64_t margin, float64_t bound, int32_t restart) 00302 { 00303 /* the knobs to be tuned ... */ 00304 /* float64_t margin = -0.95; we will go up to 95% of the 00305 distance between old variables and zero */ 00306 /* float64_t bound = 10; preset value for the start. small 00307 values give good initial 00308 feasibility but may result in slow 00309 convergence afterwards: we're too 00310 close to zero */ 00311 /* to be allocated */ 00312 float64_t *workspace; 00313 float64_t *diag_h_x; 00314 float64_t *h_y; 00315 float64_t *c_x; 00316 float64_t *c_y; 00317 float64_t *h_dot_x; 00318 float64_t *rho; 00319 float64_t *nu; 00320 float64_t *tau; 00321 float64_t *sigma; 00322 float64_t *gamma_z; 00323 float64_t *gamma_s; 00324 00325 float64_t *hat_nu; 00326 float64_t *hat_tau; 00327 00328 float64_t *delta_x; 00329 float64_t *delta_y; 00330 float64_t *delta_s; 00331 float64_t *delta_z; 00332 float64_t *delta_g; 00333 float64_t *delta_t; 00334 00335 float64_t *d; 00336 00337 /* from the header - pointers into primal and dual */ 00338 float64_t *x; 00339 float64_t *y; 00340 float64_t *g; 00341 float64_t *z; 00342 float64_t *s; 00343 float64_t *t; 00344 00345 /* auxiliary variables */ 00346 float64_t b_plus_1; 00347 float64_t c_plus_1; 00348 00349 float64_t x_h_x; 00350 float64_t primal_inf; 00351 float64_t dual_inf; 00352 00353 float64_t sigfig; 00354 float64_t primal_obj, dual_obj; 00355 float64_t mu; 00356 float64_t alfa=-1; 00357 int32_t counter = 0; 00358 00359 int32_t status = STILL_RUNNING; 00360 00361 int32_t i,j; 00362 00363 /* memory allocation */ 00364 workspace = SG_MALLOC(float64_t, (n*(m+2)+2*m)); 00365 diag_h_x = SG_MALLOC(float64_t, n); 00366 h_y = SG_MALLOC(float64_t, m*m); 00367 c_x = SG_MALLOC(float64_t, n); 00368 c_y = SG_MALLOC(float64_t, m); 00369 h_dot_x = SG_MALLOC(float64_t, n); 00370 00371 rho = SG_MALLOC(float64_t, m); 00372 nu = SG_MALLOC(float64_t, n); 00373 tau = SG_MALLOC(float64_t, n); 00374 sigma = SG_MALLOC(float64_t, n); 00375 00376 gamma_z = SG_MALLOC(float64_t, n); 00377 gamma_s = SG_MALLOC(float64_t, n); 00378 00379 hat_nu = SG_MALLOC(float64_t, n); 00380 hat_tau = SG_MALLOC(float64_t, n); 00381 00382 delta_x = SG_MALLOC(float64_t, n); 00383 delta_y = SG_MALLOC(float64_t, m); 00384 delta_s = SG_MALLOC(float64_t, n); 00385 delta_z = SG_MALLOC(float64_t, n); 00386 delta_g = SG_MALLOC(float64_t, n); 00387 delta_t = SG_MALLOC(float64_t, n); 00388 00389 d = SG_MALLOC(float64_t, n); 00390 00391 /* pointers into the external variables */ 00392 x = primal; /* n */ 00393 g = x + n; /* n */ 00394 t = g + n; /* n */ 00395 00396 y = dual; /* m */ 00397 z = y + m; /* n */ 00398 s = z + n; /* n */ 00399 00400 /* initial settings */ 00401 b_plus_1 = 1; 00402 c_plus_1 = 0; 00403 for (i=0; i<n; i++) c_plus_1 += c[i]; 00404 00405 /* get diagonal terms */ 00406 for (i=0; i<n; i++) diag_h_x[i] = h_x[(n+1)*i]; 00407 00408 /* starting point */ 00409 if (restart == 1) { 00410 /* x, y already preset */ 00411 for (i=0; i<n; i++) { /* compute g, t for primal feasibility */ 00412 g[i] = CMath::max(CMath::abs(x[i] - l[i]), bound); 00413 t[i] = CMath::max(CMath::abs(u[i] - x[i]), bound); 00414 } 00415 00416 matrix_vector(n, h_x, x, h_dot_x); /* h_dot_x = h_x * x */ 00417 00418 for (i=0; i<n; i++) { /* sigma is a dummy variable to calculate z, s */ 00419 sigma[i] = c[i] + h_dot_x[i]; 00420 for (j=0; j<m; j++) 00421 sigma[i] -= a[n*j + i] * y[j]; 00422 00423 if (sigma[i] > 0) { 00424 s[i] = bound; 00425 z[i] = sigma[i] + bound; 00426 } 00427 else { 00428 s[i] = bound - sigma[i]; 00429 z[i] = bound; 00430 } 00431 } 00432 } 00433 else { /* use default start settings */ 00434 for (i=0; i<m; i++) 00435 for (j=i; j<m; j++) 00436 h_y[i*m + j] = (i==j) ? 1 : 0; 00437 00438 for (i=0; i<n; i++) { 00439 c_x[i] = c[i]; 00440 h_x[(n+1)*i] += 1; 00441 } 00442 00443 for (i=0; i<m; i++) 00444 c_y[i] = b[i]; 00445 00446 /* and solve the system [-H_x A'; A H_y] [x, y] = [c_x; c_y] */ 00447 solve_reduced(n, m, h_x, h_y, a, x, y, c_x, c_y, workspace, 00448 PREDICTOR); 00449 00450 /* initialize the other variables */ 00451 for (i=0; i<n; i++) { 00452 g[i] = CMath::max(CMath::abs(x[i] - l[i]), bound); 00453 z[i] = CMath::max(CMath::abs(x[i]), bound); 00454 t[i] = CMath::max(CMath::abs(u[i] - x[i]), bound); 00455 s[i] = CMath::max(CMath::abs(x[i]), bound); 00456 } 00457 } 00458 00459 for (i=0, mu=0; i<n; i++) 00460 mu += z[i] * g[i] + s[i] * t[i]; 00461 mu = mu / (2*n); 00462 00463 /* the main loop */ 00464 if (verb >= STATUS) { 00465 SG_SDEBUG("counter | pri_inf | dual_inf | pri_obj | dual_obj | "); 00466 SG_SDEBUG("sigfig | alpha | nu \n"); 00467 SG_SDEBUG("-------------------------------------------------------"); 00468 SG_SDEBUG("---------------------------\n"); 00469 } 00470 00471 while (status == STILL_RUNNING) { 00472 /* predictor */ 00473 00474 /* put back original diagonal values */ 00475 for (i=0; i<n; i++) 00476 h_x[(n+1) * i] = diag_h_x[i]; 00477 00478 matrix_vector(n, h_x, x, h_dot_x); /* compute h_dot_x = h_x * x */ 00479 00480 for (i=0; i<m; i++) { 00481 rho[i] = b[i]; 00482 for (j=0; j<n; j++) 00483 rho[i] -= a[n*i + j] * x[j]; 00484 } 00485 00486 for (i=0; i<n; i++) { 00487 nu[i] = l[i] - x[i] + g[i]; 00488 tau[i] = u[i] - x[i] - t[i]; 00489 00490 sigma[i] = c[i] - z[i] + s[i] + h_dot_x[i]; 00491 for (j=0; j<m; j++) 00492 sigma[i] -= a[n*j + i] * y[j]; 00493 00494 gamma_z[i] = - z[i]; 00495 gamma_s[i] = - s[i]; 00496 } 00497 00498 /* instrumentation */ 00499 x_h_x = 0; 00500 primal_inf = 0; 00501 dual_inf = 0; 00502 00503 for (i=0; i<n; i++) { 00504 x_h_x += h_dot_x[i] * x[i]; 00505 primal_inf += CMath::sq(tau[i]); 00506 primal_inf += CMath::sq(nu[i]); 00507 dual_inf += CMath::sq(sigma[i]); 00508 } 00509 for (i=0; i<m; i++) 00510 primal_inf += CMath::sq(rho[i]); 00511 primal_inf = sqrt(primal_inf)/b_plus_1; 00512 dual_inf = sqrt(dual_inf)/c_plus_1; 00513 00514 primal_obj = 0.5 * x_h_x; 00515 dual_obj = -0.5 * x_h_x; 00516 for (i=0; i<n; i++) { 00517 primal_obj += c[i] * x[i]; 00518 dual_obj += l[i] * z[i] - u[i] * s[i]; 00519 } 00520 for (i=0; i<m; i++) 00521 dual_obj += b[i] * y[i]; 00522 00523 sigfig = log10(CMath::abs(primal_obj) + 1) - 00524 log10(CMath::abs(primal_obj - dual_obj)); 00525 sigfig = CMath::max(sigfig, 0.0); 00526 00527 /* the diagnostics - after we computed our results we will 00528 analyze them */ 00529 00530 if (counter > counter_max) status = ITERATION_LIMIT; 00531 if (sigfig > sigfig_max) status = OPTIMAL_SOLUTION; 00532 if (primal_inf > 10e100) status = PRIMAL_INFEASIBLE; 00533 if (dual_inf > 10e100) status = DUAL_INFEASIBLE; 00534 if ((primal_inf > 10e100) & (dual_inf > 10e100)) status = PRIMAL_AND_DUAL_INFEASIBLE; 00535 if (CMath::abs(primal_obj) > 10e100) status = PRIMAL_UNBOUNDED; 00536 if (CMath::abs(dual_obj) > 10e100) status = DUAL_UNBOUNDED; 00537 00538 /* write some nice routine to enforce the time limit if you 00539 _really_ want, however it's quite useless as you can compute 00540 the time from the maximum number of iterations as every 00541 iteration costs one cholesky decomposition plus a couple of 00542 backsubstitutions */ 00543 00544 /* generate report */ 00545 if ((verb >= FLOOD) | ((verb == STATUS) & (status != 0))) 00546 SG_SDEBUG("%7i | %.2e | %.2e | % .2e | % .2e | %6.3f | %.4f | %.2e\n", 00547 counter, primal_inf, dual_inf, primal_obj, dual_obj, 00548 sigfig, alfa, mu); 00549 00550 counter++; 00551 00552 if (status == 0) { /* we may keep on going, otherwise 00553 it'll cost one loop extra plus a 00554 messed up main diagonal of h_x */ 00555 /* intermediate variables (the ones with hat) */ 00556 for (i=0; i<n; i++) { 00557 hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i]; 00558 hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i]; 00559 /* diagonal terms */ 00560 d[i] = z[i] / g[i] + s[i] / t[i]; 00561 } 00562 00563 /* initialization before the cholesky solver */ 00564 for (i=0; i<n; i++) { 00565 h_x[(n+1)*i] = diag_h_x[i] + d[i]; 00566 c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - 00567 s[i] * hat_tau[i] / t[i]; 00568 } 00569 for (i=0; i<m; i++) { 00570 c_y[i] = rho[i]; 00571 for (j=i; j<m; j++) 00572 h_y[m*i + j] = 0; 00573 } 00574 00575 /* and do it */ 00576 if (!solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace, PREDICTOR)) 00577 { 00578 status=INCONSISTENT; 00579 goto exit_optimizer; 00580 } 00581 00582 for (i=0; i<n; i++) { 00583 /* backsubstitution */ 00584 delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i]; 00585 delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i]; 00586 00587 delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i]; 00588 delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i]; 00589 00590 /* central path (corrector) */ 00591 gamma_z[i] = mu / g[i] - z[i] - delta_z[i] * delta_g[i] / g[i]; 00592 gamma_s[i] = mu / t[i] - s[i] - delta_s[i] * delta_t[i] / t[i]; 00593 00594 /* (some more intermediate variables) the hat variables */ 00595 hat_nu[i] = nu[i] + g[i] * gamma_z[i] / z[i]; 00596 hat_tau[i] = tau[i] - t[i] * gamma_s[i] / s[i]; 00597 00598 /* initialization before the cholesky */ 00599 c_x[i] = sigma[i] - z[i] * hat_nu[i] / g[i] - s[i] * hat_tau[i] / t[i]; 00600 } 00601 00602 for (i=0; i<m; i++) { /* comput c_y and rho */ 00603 c_y[i] = rho[i]; 00604 for (j=i; j<m; j++) 00605 h_y[m*i + j] = 0; 00606 } 00607 00608 /* and do it */ 00609 solve_reduced(n, m, h_x, h_y, a, delta_x, delta_y, c_x, c_y, workspace, 00610 CORRECTOR); 00611 00612 for (i=0; i<n; i++) { 00613 /* backsubstitution */ 00614 delta_s[i] = s[i] * (delta_x[i] - hat_tau[i]) / t[i]; 00615 delta_z[i] = z[i] * (hat_nu[i] - delta_x[i]) / g[i]; 00616 00617 delta_g[i] = g[i] * (gamma_z[i] - delta_z[i]) / z[i]; 00618 delta_t[i] = t[i] * (gamma_s[i] - delta_s[i]) / s[i]; 00619 } 00620 00621 alfa = -1; 00622 for (i=0; i<n; i++) { 00623 alfa = CMath::min(alfa, delta_g[i]/g[i]); 00624 alfa = CMath::min(alfa, delta_t[i]/t[i]); 00625 alfa = CMath::min(alfa, delta_s[i]/s[i]); 00626 alfa = CMath::min(alfa, delta_z[i]/z[i]); 00627 } 00628 alfa = (margin - 1) / alfa; 00629 00630 /* compute mu */ 00631 for (i=0, mu=0; i<n; i++) 00632 mu += z[i] * g[i] + s[i] * t[i]; 00633 mu = mu / (2*n); 00634 mu = mu * CMath::sq((alfa - 1) / (alfa + 10)); 00635 00636 for (i=0; i<n; i++) { 00637 x[i] += alfa * delta_x[i]; 00638 g[i] += alfa * delta_g[i]; 00639 t[i] += alfa * delta_t[i]; 00640 z[i] += alfa * delta_z[i]; 00641 s[i] += alfa * delta_s[i]; 00642 } 00643 00644 for (i=0; i<m; i++) 00645 y[i] += alfa * delta_y[i]; 00646 } 00647 } 00648 00649 exit_optimizer: 00650 if ((status == 1) && (verb >= STATUS)) { 00651 SG_SDEBUG("----------------------------------------------------------------------------------\n"); 00652 SG_SDEBUG("optimization converged\n"); 00653 } 00654 00655 /* free memory */ 00656 SG_FREE(workspace); 00657 SG_FREE(diag_h_x); 00658 SG_FREE(h_y); 00659 SG_FREE(c_x); 00660 SG_FREE(c_y); 00661 SG_FREE(h_dot_x); 00662 00663 SG_FREE(rho); 00664 SG_FREE(nu); 00665 SG_FREE(tau); 00666 SG_FREE(sigma); 00667 SG_FREE(gamma_z); 00668 SG_FREE(gamma_s); 00669 00670 SG_FREE(hat_nu); 00671 SG_FREE(hat_tau); 00672 00673 SG_FREE(delta_x); 00674 SG_FREE(delta_y); 00675 SG_FREE(delta_s); 00676 SG_FREE(delta_z); 00677 SG_FREE(delta_g); 00678 SG_FREE(delta_t); 00679 00680 SG_FREE(d); 00681 00682 /* and return to sender */ 00683 return status; 00684 } 00685 }