SHOGUN
v1.1.0
|
00001 /****************************************************************************** 00002 *** GPDT - Gradient Projection Decomposition Technique *** 00003 ****************************************************************************** 00004 *** *** 00005 *** GPDT is a C++ software designed to train large-scale Support Vector *** 00006 *** Machines for binary classification in both scalar and distributed *** 00007 *** memory parallel environments. It uses the Joachims' problem *** 00008 *** decomposition technique to split the whole quadratic programming (QP) *** 00009 *** problem into a sequence of smaller QP subproblems, each one being *** 00010 *** solved by a suitable gradient projection method (GPM). The presently *** 00011 *** implemented GPMs are the Generalized Variable Projection Method *** 00012 *** GVPM (T. Serafini, G. Zanghirati, L. Zanni, "Gradient Projection *** 00013 *** Methods for Quadratic Programs and Applications in Training Support *** 00014 *** Vector Machines"; Optim. Meth. Soft. 20, 2005, 353-378) and the *** 00015 *** Dai-Fletcher Method DFGPM (Y. Dai and R. Fletcher,"New Algorithms for *** 00016 *** Singly Linear Constrained Quadratic Programs Subject to Lower and *** 00017 *** Upper Bounds"; Math. Prog. to appear). *** 00018 *** *** 00019 *** Authors: *** 00020 *** Thomas Serafini, Luca Zanni *** 00021 *** Dept. of Mathematics, University of Modena and Reggio Emilia - ITALY *** 00022 *** serafini.thomas@unimo.it, zanni.luca@unimo.it *** 00023 *** Gaetano Zanghirati *** 00024 *** Dept. of Mathematics, University of Ferrara - ITALY *** 00025 *** g.zanghirati@unife.it *** 00026 *** *** 00027 *** Software homepage: http://dm.unife.it/gpdt *** 00028 *** *** 00029 *** This work is supported by the Italian FIRB Projects *** 00030 *** 'Statistical Learning: Theory, Algorithms and Applications' *** 00031 *** (grant RBAU01877P), http://slipguru.disi.unige.it/ASTA *** 00032 *** and *** 00033 *** 'Parallel Algorithms and Numerical Nonlinear Optimization' *** 00034 *** (grant RBAU01JYPN), http://dm.unife.it/pn2o *** 00035 *** *** 00036 *** Copyright (C) 2004-2008 by T. Serafini, G. Zanghirati, L. Zanni. *** 00037 *** *** 00038 *** COPYRIGHT NOTIFICATION *** 00039 *** *** 00040 *** Permission to copy and modify this software and its documentation *** 00041 *** for internal research use is granted, provided that this notice is *** 00042 *** retained thereon and on all copies or modifications. The authors and *** 00043 *** their respective Universities makes no representations as to the *** 00044 *** suitability and operability of this software for any purpose. It is *** 00045 *** provided "as is" without express or implied warranty. *** 00046 *** Use of this software for commercial purposes is expressly prohibited *** 00047 *** without contacting the authors. *** 00048 *** *** 00049 *** This program is free software; you can redistribute it and/or modify *** 00050 *** it under the terms of the GNU General Public License as published by *** 00051 *** the Free Software Foundation; either version 3 of the License, or *** 00052 *** (at your option) any later version. *** 00053 *** *** 00054 *** This program is distributed in the hope that it will be useful, *** 00055 *** but WITHOUT ANY WARRANTY; without even the implied warranty of *** 00056 *** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *** 00057 *** GNU General Public License for more details. *** 00058 *** *** 00059 *** You should have received a copy of the GNU General Public License *** 00060 *** along with this program; if not, write to the Free Software *** 00061 *** Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. *** 00062 *** *** 00063 *** File: gpdtsolve.cpp *** 00064 *** Type: scalar *** 00065 *** Version: 1.0 *** 00066 *** Date: November, 2006 *** 00067 *** Revision: 2 *** 00068 *** *** 00069 *** SHOGUN adaptions Written (W) 2006-2009 Soeren Sonnenburg *** 00070 ******************************************************************************/ 00071 #include <math.h> 00072 #include <stdio.h> 00073 #include <stdlib.h> 00074 #include <string.h> 00075 #include <time.h> 00076 00077 #include <shogun/classifier/svm/gpm.h> 00078 #include <shogun/classifier/svm/gpdt.h> 00079 #include <shogun/classifier/svm/gpdtsolve.h> 00080 #include <shogun/lib/Signal.h> 00081 #include <shogun/io/SGIO.h> 00082 00083 using namespace shogun; 00084 00085 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00086 namespace shogun 00087 { 00088 #define y_in(i) y[index_in[(i)]] 00089 #define y_out(i) y[index_out[(i)]] 00090 #define alpha_in(i) alpha[index_in[(i)]] 00091 #define alpha_out(i) alpha[index_out[(i)]] 00092 #define minfty (-1.0e+30) // minus infinity 00093 00094 uint32_t Randnext = 1; 00095 00096 #define ThRand (Randnext = Randnext * 1103515245L + 12345L) 00097 #define ThRandPos ((Randnext = Randnext * 1103515245L + 12345L) & 0x7fffffff) 00098 00099 FILE *fp; 00100 00101 /* utility routines prototyping */ 00102 void quick_si (int32_t a[], int32_t k); 00103 void quick_s3 (int32_t a[], int32_t k, int32_t ia[]); 00104 void quick_s2 (float64_t a[], int32_t k, int32_t ia[]); 00105 00106 /******************************************************************************/ 00107 /*** Class for caching strategy implementation ***/ 00108 /******************************************************************************/ 00109 class sCache 00110 { 00111 00112 public: 00113 sCache (sKernel* sk, int32_t Mbyte, int32_t ell); 00114 ~sCache (); 00115 00116 cachetype *FillRow (int32_t row, int32_t IsC = 0); 00117 cachetype *GetRow (int32_t row); 00118 00119 int32_t DivideMP (int32_t *out, int32_t *in, int32_t n); 00120 00121 /*** Itarations counter ***/ 00122 void Iteration() { nit++; } 00123 00124 /*** Cache size control ***/ 00125 int32_t CheckCycle() 00126 { 00127 int32_t us; 00128 cache_entry *pt = first_free->next; 00129 00130 for (us = 0; pt != first_free; us++, pt = pt->next); 00131 if (us != maxmw-1) 00132 return 1; 00133 else 00134 return 0; 00135 } 00136 00137 private: 00138 00139 struct cache_entry 00140 { 00141 int32_t row; // unused row 00142 int32_t last_access_it; 00143 cache_entry *prev, *next; 00144 cachetype *data; 00145 }; 00146 00147 sKernel* KER; 00148 int32_t maxmw, ell; 00149 int32_t nit; 00150 00151 cache_entry *mw; 00152 cache_entry *first_free; 00153 cache_entry **pindmw; // 0 if unused row 00154 cachetype *onerow; 00155 00156 cachetype *FindFree(int32_t row, int32_t IsC); 00157 }; 00158 00159 00160 /******************************************************************************/ 00161 /*** Cache class constructor ***/ 00162 /******************************************************************************/ 00163 sCache::sCache(sKernel* sk, int32_t Mbyte, int32_t _ell) : KER(sk), ell(_ell) 00164 { 00165 int32_t i; 00166 00167 // size in dwords of one cache row 00168 maxmw = (sizeof(cache_entry) + sizeof(cache_entry *) 00169 + ell*sizeof(cachetype)) / 4; 00170 // number of cache rows 00171 maxmw = Mbyte*1024*(1024/4) / maxmw; 00172 00173 /* arrays allocation */ 00174 mw = SG_MALLOC(cache_entry, maxmw); 00175 pindmw = SG_MALLOC(cache_entry*, ell); 00176 onerow = SG_MALLOC(cachetype, ell); 00177 00178 /* arrays initialization */ 00179 for (i = 0; i < maxmw; i++) 00180 { 00181 mw[i].prev = (i == 0 ? &mw[maxmw-1] : &mw[i-1]); 00182 mw[i].next = (i == maxmw-1 ? &mw[0] : &mw[i+1]); 00183 mw[i].data = SG_MALLOC(cachetype, ell); 00184 mw[i].row = -1; // unused row 00185 mw[i].last_access_it = -1; 00186 } 00187 for (i = 0; i < ell; i++) 00188 pindmw[i] = 0; 00189 00190 first_free = &mw[0]; 00191 nit = 0; 00192 } 00193 00194 /******************************************************************************/ 00195 /*** Cache class destructor ***/ 00196 /******************************************************************************/ 00197 sCache::~sCache() 00198 { 00199 int32_t i; 00200 00201 for (i = maxmw-1; i >= 0; i--) 00202 SG_FREE(mw[i].data); 00203 00204 SG_FREE(onerow); 00205 SG_FREE(pindmw); 00206 SG_FREE(mw); 00207 } 00208 00209 00210 /******************************************************************************/ 00211 /*** Retrieve a cached row ***/ 00212 /******************************************************************************/ 00213 cachetype *sCache::GetRow(int32_t row) 00214 { 00215 cache_entry *c; 00216 00217 c = pindmw[row]; 00218 if (c == NULL) 00219 return NULL; 00220 00221 c->last_access_it = nit; 00222 if (c == first_free) 00223 { 00224 first_free = first_free->next; 00225 } 00226 else 00227 { 00228 // move "row" as the most recently used. 00229 c->prev->next = c->next; 00230 c->next->prev = c->prev; 00231 c->next = first_free; 00232 c->prev = first_free->prev; 00233 first_free->prev = c; 00234 c->prev->next = c; 00235 } 00236 00237 return c->data; 00238 } 00239 00240 /****************************************************************************** 00241 *** Look for a free cache row *** 00242 *** IMPORTANT: call this method only if you are sure that "row" *** 00243 *** is not already in the cache ( i.e. after calling GetRow() ) *** 00244 ******************************************************************************/ 00245 cachetype *sCache::FindFree(int32_t row, int32_t IsC) 00246 { 00247 cachetype *pt; 00248 00249 if (first_free->row != -1) // cache row already contains data 00250 { 00251 if (first_free->last_access_it == nit || IsC) 00252 return 0; 00253 else 00254 pindmw[first_free->row] = 0; 00255 } 00256 first_free->row = row; 00257 first_free->last_access_it = nit; 00258 pindmw[row] = first_free; 00259 00260 pt = first_free->data; 00261 first_free = first_free->next; 00262 00263 return pt; 00264 } 00265 00266 00267 /******************************************************************************/ 00268 /*** Enter data in a cache row ***/ 00269 /******************************************************************************/ 00270 cachetype *sCache::FillRow(int32_t row, int32_t IsC) 00271 { 00272 int32_t j; 00273 cachetype *pt; 00274 00275 pt = GetRow(row); 00276 if (pt != NULL) 00277 return pt; 00278 00279 pt = FindFree(row, IsC); 00280 if (pt == 0) 00281 pt = onerow; 00282 00283 // Compute all the row elements 00284 for (j = 0; j < ell; j++) 00285 pt[j] = (cachetype)KER->Get(row, j); 00286 return pt; 00287 } 00288 00289 00290 /******************************************************************************/ 00291 /*** Expand a sparse row in a full cache row ***/ 00292 /******************************************************************************/ 00293 int32_t sCache::DivideMP(int32_t *out, int32_t *in, int32_t n) 00294 { 00295 /******************************************************************** 00296 * Input meaning: * 00297 * in = vector containing row to be extracted in the cache * 00298 * n = size of in * 00299 * out = the indexes of "in" of the components to be computed * 00300 * by this processor (first those in the cache, then the * 00301 * ones not yet computed) * 00302 * Returns: the number of components of this processor * 00303 ********************************************************************/ 00304 00305 int32_t *remained, nremained, k; 00306 int32_t i; 00307 00308 remained = SG_MALLOC(int32_t, n); 00309 00310 nremained = 0; 00311 k = 0; 00312 for (i = 0; i < n; i++) 00313 { 00314 if (pindmw[in[i]] != NULL) 00315 out[k++] = i; 00316 else 00317 remained[nremained++] = i; 00318 } 00319 for (i = 0; i < nremained; i++) 00320 out[k++] = remained[i]; 00321 00322 SG_FREE(remained); 00323 return n; 00324 } 00325 00326 /******************************************************************************/ 00327 /*** Check solution optimality ***/ 00328 /******************************************************************************/ 00329 int32_t QPproblem::optimal() 00330 { 00331 /*********************************************************************** 00332 * Returns 1 if the computed solution is optimal, otherwise returns 0. * 00333 * To verify the optimality it checks the KKT optimality conditions. * 00334 ***********************************************************************/ 00335 register int32_t i, j, margin_sv_number, z, k, s, kin, z1, znew=0, nnew; 00336 00337 float64_t gx_i, aux, s1, s2; 00338 00339 /* sort -y*grad and store in ing the swaps done */ 00340 for (j = 0; j < ell; j++) 00341 { 00342 grad[j] = y[j] - st[j]; 00343 ing[j] = j; 00344 } 00345 00346 quick_s2(grad,ell,ing); 00347 00348 /* compute bee */ 00349 margin_sv_number = 0; 00350 00351 for (i = chunk_size - 1; i >= 0; i--) 00352 if (is_free(index_in[i])) 00353 { 00354 margin_sv_number++; 00355 bee = y_in(i) - st[index_in[i]]; 00356 break; 00357 } 00358 00359 if (margin_sv_number > 0) 00360 { 00361 aux=-1.0; 00362 for (i = nb-1; i >= 0; i--) 00363 { 00364 gx_i = bee + st[index_out[i]]; 00365 if ((is_zero(index_out[i]) && (gx_i*y_out(i) < (1.0-delta))) || 00366 (is_bound(index_out[i]) && (gx_i*y_out(i) > (1.0+delta))) || 00367 (is_free(index_out[i]) && 00368 ((gx_i*y_out(i) < 1.0-delta) || (gx_i*y_out(i) > 1.0+delta)))) 00369 { 00370 if (fabs(gx_i*y_out(i) - 1.0) > aux) 00371 aux = fabs(gx_i*y_out(i) - 1.0); 00372 } 00373 } 00374 } 00375 else 00376 { 00377 for (i = ell - 1; i >= 0; i--) 00378 if (is_free(i)) 00379 { 00380 margin_sv_number++; 00381 bee = y[i] - st[i]; 00382 break; 00383 } 00384 if (margin_sv_number == 0) 00385 { 00386 s1 = -minfty; 00387 s2 = -minfty; 00388 for (j = 0; j < ell; j++) 00389 if ( (alpha[ing[j]] > DELTAsv) && (y[ing[j]] >= 0) ) 00390 { 00391 s1 = grad[j]; 00392 break; 00393 } 00394 for (j = 0; j < ell; j++) 00395 if ( (alpha[ing[j]] < c_const-DELTAsv) && (y[ing[j]] <= 0) ) 00396 { 00397 s2 = grad[j]; 00398 break; 00399 } 00400 if (s1 < s2) 00401 aux = s1; 00402 else 00403 aux = s2; 00404 00405 s1 = minfty; 00406 s2 = minfty; 00407 for (j = ell-1; j >=0; j--) 00408 if ( (alpha[ing[j]] > DELTAsv) && (y[ing[j]] <= 0) ) 00409 { 00410 s1 = grad[j]; 00411 break; 00412 } 00413 for (j = ell-1; j >=0; j--) 00414 if ( (alpha[ing[j]] < c_const-DELTAsv) && (y[ing[j]] >= 0) ) 00415 { 00416 s2 = grad[j]; 00417 break; 00418 } 00419 if (s2 > s1) s1 = s2; 00420 00421 bee = 0.5 * (s1+aux); 00422 } 00423 00424 aux = -1.0; 00425 for (i = ell-1; i >= 0; i--) 00426 { 00427 gx_i = bee + st[i]; 00428 if ((is_zero(i) && (gx_i*y[i] < (1.0-delta))) || 00429 (is_bound(i) && (gx_i*y[i] > (1.0+delta))) || 00430 (is_free(i) && 00431 ((gx_i*y[i] < 1.0-delta) || (gx_i*y[i] > 1.0+delta)))) 00432 { 00433 if (fabs(gx_i*y[i] - 1.0) > aux) 00434 aux = fabs(gx_i*y[i] - 1.0); 00435 } 00436 } 00437 } 00438 00439 if (aux < 0.0) 00440 return 1; 00441 else 00442 { 00443 if (verbosity > 1) 00444 SG_INFO(" Max KKT violation: %lf\n", aux); 00445 else if (verbosity > 0) 00446 SG_INFO(" %lf\n", aux); 00447 00448 if (fabs(kktold-aux) < delta*0.01 && aux < delta*2.0) 00449 { 00450 if (DELTAvpm > InitialDELTAvpm*0.1) 00451 { 00452 DELTAvpm = (DELTAvpm*0.5 > InitialDELTAvpm*0.1 ? 00453 DELTAvpm*0.5 : InitialDELTAvpm*0.1); 00454 SG_INFO("Inner tolerance changed to: %lf\n", DELTAvpm); 00455 } 00456 } 00457 00458 kktold = aux; 00459 00460 /***************************************************************************** 00461 *** Update the working set (T. Serafini, L. Zanni, "On the Working Set *** 00462 *** Selection in Gradient Projection-based Decomposition Techniques for *** 00463 *** Support Vector Machines"; Optim. Meth. Soft. 20, 2005). *** 00464 *****************************************************************************/ 00465 for (j = 0; j < chunk_size; j++) 00466 inold[j] = index_in[j]; 00467 00468 z = -1; /* index of the last entered component from the top */ 00469 z1 = ell; /* index of the last entered component from the bottom */ 00470 k = 0; 00471 j = 0; 00472 while (k < q) 00473 { 00474 i = z + 1; /* index of the candidate components from the top */ 00475 while (i < z1) 00476 { 00477 if ( is_free(ing[i]) || 00478 (-y[ing[i]]>=0 && is_zero(ing[i])) || 00479 (-y[ing[i]]<=0 && is_bound(ing[i])) 00480 ) 00481 { 00482 znew = i; /* index of the component to select from the top */ 00483 break; 00484 } 00485 i++; 00486 } 00487 if (i == z1) break; 00488 00489 s = z1 - 1; 00490 while (znew < s) 00491 { 00492 if ( is_free(ing[s]) || 00493 (y[ing[s]]>=0 && is_zero(ing[s])) || 00494 (y[ing[s]]<=0 && is_bound(ing[s])) 00495 ) 00496 { 00497 z1 = s; 00498 z = znew; 00499 break; 00500 } 00501 s--; 00502 } 00503 if (znew == s) break; 00504 00505 index_in[k++] = ing[z]; 00506 index_in[k++] = ing[z1]; 00507 } 00508 00509 if (k < q) 00510 { 00511 if (verbosity > 1) 00512 SG_INFO(" New q: %i\n", k); 00513 q = k; 00514 } 00515 00516 quick_si(index_in, q); 00517 00518 s = 0; 00519 j = 0; 00520 for (k = 0; k < chunk_size; k++) 00521 { 00522 z = inold[k]; 00523 for (i = j; i < q; i++) 00524 if (z <= index_in[i]) 00525 break; 00526 00527 if (i == q) 00528 { 00529 for (i = k; i < chunk_size; i++) 00530 { 00531 ing[s] = inold[i]; /* older components not in the new basis */ 00532 s = s+1; 00533 } 00534 break; 00535 } 00536 00537 if (z == index_in[i]) 00538 j = i + 1; /* the component is still in the basis */ 00539 else 00540 { 00541 ing[s] = z; /* older component not in the new basis */ 00542 s = s + 1; 00543 j = i; 00544 } 00545 } 00546 00547 for (i = 0; i < s; i++) 00548 { 00549 bmemrid[i] = bmem[ing[i]]; 00550 pbmr[i] = i; 00551 } 00552 00553 quick_s3(bmemrid, s, pbmr); 00554 00555 /* check if support vectors not at bound enter the basis */ 00556 j = q; 00557 i = 0; 00558 while (j < chunk_size && i < s) 00559 { 00560 if (is_free(ing[pbmr[i]])) 00561 { 00562 index_in[j] = ing[pbmr[i]]; 00563 j++; 00564 } 00565 i++; 00566 } 00567 00568 /* choose which bound variables keep in basis (alpha = 0 or alpha = C) */ 00569 if (j < chunk_size) 00570 { 00571 i = 0; 00572 while (j < chunk_size && i < s) 00573 { 00574 if (is_zero(ing[pbmr[i]])) 00575 { 00576 index_in[j] = ing[pbmr[i]]; 00577 j++; 00578 } 00579 i++; 00580 } 00581 if (j < chunk_size) 00582 { 00583 i = 0; 00584 while (j < chunk_size && i < s) 00585 { 00586 if (is_bound(ing[pbmr[i]])) 00587 { 00588 index_in[j] = ing[pbmr[i]]; 00589 j++; 00590 } 00591 i++; 00592 } 00593 } 00594 } 00595 00596 quick_si(index_in, chunk_size); 00597 00598 for (i = 0; i < chunk_size; i++) 00599 bmem[index_in[i]]++; 00600 00601 j = 0; 00602 k = 0; 00603 for (i = 0; i < chunk_size; i++) 00604 { 00605 for (z = j; z < index_in[i]; z++) 00606 { 00607 index_out[k] = z; 00608 k++; 00609 } 00610 j = index_in[i]+1; 00611 } 00612 for (z = j; z < ell; z++) 00613 { 00614 index_out[k] = z; 00615 k++; 00616 } 00617 00618 for (i = 0; i < nb; i++) 00619 bmem[index_out[i]] = 0; 00620 00621 kin = 0; 00622 j = 0; 00623 for (k = 0; k < chunk_size; k++) 00624 { 00625 z = index_in[k]; 00626 for (i = j; i < chunk_size; i++) 00627 if (z <= inold[i]) 00628 break; 00629 if (i == chunk_size) 00630 { 00631 for (s = k; s < chunk_size; s++) 00632 { 00633 incom[s] = -1; 00634 cec[index_in[s]]++; 00635 } 00636 kin = kin + chunk_size - k ; 00637 break; 00638 } 00639 00640 if (z == inold[i]) 00641 { 00642 incom[k] = i; 00643 j = i+1; 00644 } 00645 else 00646 { 00647 incom[k] = -1; 00648 j = i; 00649 kin = kin + 1; 00650 cec[index_in[k]]++; 00651 } 00652 } 00653 00654 nnew = kin & (~1); 00655 if (nnew < 10) 00656 nnew = 10; 00657 if (nnew < chunk_size/10) 00658 nnew = chunk_size/10; 00659 if (nnew < q) 00660 { 00661 q = nnew; 00662 q = q & (~1); 00663 } 00664 00665 if (kin == 0) 00666 { 00667 DELTAkin *= 0.1; 00668 if (DELTAkin < 1.0e-6) 00669 { 00670 SG_INFO("\n***ERROR***: GPDT stops with tolerance"); 00671 SG_INFO( 00672 " %lf because it is unable to change the working set.\n", kktold); 00673 return 1; 00674 } 00675 else 00676 { 00677 SG_INFO("Inner tolerance temporary changed to:"); 00678 SG_INFO(" %e\n", DELTAvpm*DELTAkin); 00679 } 00680 } 00681 else 00682 DELTAkin = 1.0; 00683 00684 if (verbosity > 1) 00685 { 00686 SG_INFO(" Working set: new components: %i", kin); 00687 SG_INFO(", new parameter n: %i\n", q); 00688 } 00689 00690 return 0; 00691 } 00692 } 00693 00694 /******************************************************************************/ 00695 /*** Optional preprocessing: random distribution ***/ 00696 /******************************************************************************/ 00697 int32_t QPproblem::Preprocess0(int32_t *aux, int32_t *sv) 00698 { 00699 int32_t i, j; 00700 00701 Randnext = 1; 00702 memset(sv, 0, ell*sizeof(int32_t)); 00703 for (i = 0; i < chunk_size; i++) 00704 { 00705 do 00706 { 00707 j = ThRandPos % ell; 00708 } while (sv[j] != 0); 00709 sv[j] = 1; 00710 } 00711 return(chunk_size); 00712 } 00713 00714 /******************************************************************************/ 00715 /*** Optional preprocessing: block parallel distribution ***/ 00716 /******************************************************************************/ 00717 int32_t QPproblem::Preprocess1(sKernel* kernel, int32_t *aux, int32_t *sv) 00718 { 00719 int32_t s; // elements owned by the processor 00720 int32_t sl; // elements of the n-1 subproblems 00721 int32_t n, i, off, j, k, ll; 00722 int32_t nsv, nbsv; 00723 int32_t *sv_loc, *bsv_loc, *sp_y; 00724 float32_t *sp_D=NULL; 00725 float64_t *sp_alpha, *sp_h; 00726 00727 s = ell; 00728 /* divide the s elements into n blocks smaller than preprocess_size */ 00729 n = (s + preprocess_size - 1) / preprocess_size; 00730 sl = 1 + s / n; 00731 00732 if (verbosity > 0) 00733 { 00734 SG_INFO(" Preprocessing: examples = %d", s); 00735 SG_INFO(", subp. = %d", n); 00736 SG_INFO(", size = %d\n",sl); 00737 } 00738 00739 sv_loc = SG_MALLOC(int32_t, s); 00740 bsv_loc = SG_MALLOC(int32_t, s); 00741 sp_alpha = SG_MALLOC(float64_t, sl); 00742 sp_h = SG_MALLOC(float64_t, sl); 00743 sp_y = SG_MALLOC(int32_t, sl); 00744 00745 if (sl < 500) 00746 sp_D = SG_MALLOC(float32_t, sl*sl); 00747 00748 for (i = 0; i < sl; i++) 00749 sp_h[i] = -1.0; 00750 memset(alpha, 0, ell*sizeof(float64_t)); 00751 00752 /* randomly reorder the component to select */ 00753 for (i = 0; i < ell; i++) 00754 aux[i] = i; 00755 Randnext = 1; 00756 for (i = 0; i < ell; i++) 00757 { 00758 j = ThRandPos % ell; 00759 k = ThRandPos % ell; 00760 ll = aux[j]; aux[j] = aux[k]; aux[k] = ll; 00761 } 00762 00763 nbsv = nsv = 0; 00764 for (i = 0; i < n; i++) 00765 { 00766 if (verbosity > 0) 00767 SG_INFO("%d...", i); 00768 SplitParts(s, i, n, &ll, &off); 00769 00770 if (sl < 500) 00771 { 00772 for (j = 0; j < ll; j++) 00773 { 00774 sp_y[j] = y[aux[j+off]]; 00775 for (k = j; k < ll; k++) 00776 sp_D[k*sl + j] = sp_D[j*sl + k] 00777 = y[aux[j+off]] * y[aux[k+off]] 00778 * (float32_t)kernel->Get(aux[j+off], aux[k+off]); 00779 } 00780 00781 memset(sp_alpha, 0, sl*sizeof(float64_t)); 00782 00783 /* call the gradient projection QP solver */ 00784 gpm_solver(projection_solver, projection_projector, ll, sp_D, sp_h, 00785 c_const, 0.0, sp_y, sp_alpha, delta*10, NULL); 00786 } 00787 else 00788 { 00789 QPproblem p2; 00790 p2.Subproblem(*this, ll, aux + off); 00791 p2.chunk_size = (int32_t) ((float64_t)chunk_size / sqrt((float64_t)n)); 00792 p2.q = (int32_t) ((float64_t)q / sqrt((float64_t)n)); 00793 p2.maxmw = ll*ll*4 / (1024 * 1024); 00794 if (p2.maxmw > maxmw / 2) 00795 p2.maxmw = maxmw / 2; 00796 p2.verbosity = 0; 00797 p2.delta = delta * 10.0; 00798 p2.PreprocessMode = 0; 00799 kernel->KernelEvaluations += p2.gpdtsolve(sp_alpha); 00800 } 00801 00802 for (j = 0; j < ll; j++) 00803 { 00804 /* modify bound support vector approximation */ 00805 if (sp_alpha[j] < (c_const-DELTAsv)) 00806 sp_alpha[j] = 0.0; 00807 else 00808 sp_alpha[j] = c_const; 00809 if (sp_alpha[j] > DELTAsv) 00810 { 00811 if (sp_alpha[j] < (c_const-DELTAsv)) 00812 sv_loc[nsv++] = aux[j+off]; 00813 else 00814 bsv_loc[nbsv++] = aux[j+off]; 00815 alpha[aux[j+off]] = sp_alpha[j]; 00816 } 00817 } 00818 } 00819 00820 Randnext = 1; 00821 00822 /* add the known support vectors to the working set */ 00823 memset(sv, 0, ell*sizeof(int32_t)); 00824 ll = (nsv < chunk_size ? nsv : chunk_size); 00825 for (i = 0; i < ll; i++) 00826 { 00827 do { 00828 j = sv_loc[ThRandPos % nsv]; 00829 } while (sv[j] != 0); 00830 sv[j] = 1; 00831 } 00832 00833 /* add the known bound support vectors to the working set */ 00834 ll = ((nsv+nbsv) < chunk_size ? (nsv+nbsv) : chunk_size); 00835 for (; i < ll; i++) 00836 { 00837 do { 00838 j = bsv_loc[ThRandPos % nbsv]; 00839 } while (sv[j] != 0); 00840 sv[j] = 1; 00841 } 00842 00843 /* eventually fill up the working set with other components 00844 randomly chosen */ 00845 for (; i < chunk_size; i++) 00846 { 00847 do { 00848 j = ThRandPos % ell; 00849 } while (sv[j] != 0); 00850 sv[j] = 1; 00851 } 00852 00853 00854 /* dealloc temporary arrays */ 00855 if (sl < 500) SG_FREE(sp_D); 00856 SG_FREE(sp_y ); 00857 SG_FREE(sp_h ); 00858 SG_FREE(sv_loc ); 00859 SG_FREE(bsv_loc ); 00860 SG_FREE(sp_alpha); 00861 00862 if (verbosity > 0) 00863 { 00864 SG_INFO("\n Preprocessing: SV = %d", nsv); 00865 SG_INFO(", BSV = %d\n", nbsv); 00866 } 00867 00868 return(nsv); 00869 } 00870 00871 /******************************************************************************/ 00872 /*** Compute the QP problem solution ***/ 00873 /******************************************************************************/ 00874 float64_t QPproblem::gpdtsolve(float64_t *solution) 00875 { 00876 int32_t i, j, k, z, iin, jin, nit, tot_vpm_iter, lsCount; 00877 int32_t tot_vpm_secant, projCount, proximal_count; 00878 int32_t vpmWarningThreshold; 00879 int32_t nzin, nzout; 00880 int32_t *sp_y; /* labels vector */ 00881 int32_t *indnzin, *indnzout; /* nonzero components indices vectors */ 00882 float32_t *sp_D; /* quadratic part of the objective function */ 00883 float64_t *sp_h, *sp_hloc, /* linear part of the objective function */ 00884 *sp_alpha,*stloc; /* variables and gradient updating vectors */ 00885 float64_t sp_e, aux, fval, tau_proximal_this, dfval; 00886 float64_t *vau; 00887 float64_t *weight; 00888 float64_t tot_prep_time, tot_vpm_time, tot_st_time, total_time; 00889 sCache *Cache; 00890 cachetype *ptmw; 00891 clock_t t, ti; 00892 00893 Cache = new sCache(KER, maxmw, ell); 00894 if (chunk_size > ell) chunk_size = ell; 00895 00896 if (chunk_size <= 20) 00897 vpmWarningThreshold = 30*chunk_size; 00898 else if (chunk_size <= 200) 00899 vpmWarningThreshold = 20*chunk_size + 200; 00900 else 00901 vpmWarningThreshold = 10*chunk_size + 2200; 00902 00903 kktold = 10000.0; 00904 if (delta <= 5e-3) 00905 { 00906 if ( (chunk_size <= 20) | ((float64_t)chunk_size/ell <= 0.001) ) 00907 DELTAvpm = delta * 0.1; 00908 else if ( (chunk_size <= 200) | ((float64_t)chunk_size/ell <= 0.005) ) 00909 DELTAvpm = delta * 0.5; 00910 else 00911 DELTAvpm = delta; 00912 } 00913 else 00914 { 00915 if ( (chunk_size <= 200) | ((float64_t)chunk_size/ell <= 0.005) ) 00916 DELTAvpm = (1e-3 < delta*0.1) ? 1e-3 : delta*0.1; 00917 else 00918 DELTAvpm = 5e-3; 00919 } 00920 00921 InitialDELTAvpm = DELTAvpm; 00922 DELTAsv = EPS_SV * c_const; 00923 DELTAkin = 1.0; 00924 00925 q = q & (~1); 00926 nb = ell - chunk_size; 00927 tot_vpm_iter = 0; 00928 tot_vpm_secant = 0; 00929 00930 tot_prep_time = tot_vpm_time = tot_st_time = total_time = 0.0; 00931 00932 ti = clock(); 00933 00934 /* arrays allocation */ 00935 SG_DEBUG("ell:%d, chunk_size:%d, nb:%d dim:%d\n", ell, chunk_size,nb, dim); 00936 ing = SG_MALLOC(int32_t, ell); 00937 inaux = SG_MALLOC(int32_t, ell); 00938 index_in = SG_MALLOC(int32_t, chunk_size); 00939 index_out = SG_MALLOC(int32_t, ell); 00940 indnzout = SG_MALLOC(int32_t, nb); 00941 alpha = SG_MALLOC(float64_t, ell); 00942 00943 memset(alpha, 0, ell*sizeof(float64_t)); 00944 memset(ing, 0, ell*sizeof(int32_t)); 00945 00946 if (verbosity > 0 && PreprocessMode != 0) 00947 SG_INFO("\n*********** Begin setup step...\n"); 00948 t = clock(); 00949 00950 switch(PreprocessMode) 00951 { 00952 case 1: Preprocess1(KER, inaux, ing); break; 00953 case 0: 00954 default: 00955 Preprocess0(inaux, ing); break; 00956 } 00957 00958 for (j = k = i = 0; i < ell; i++) 00959 if (ing[i] == 0) 00960 index_out[j++] = i; 00961 else 00962 index_in[k++] = i; 00963 00964 t = clock() - t; 00965 if (verbosity > 0 && PreprocessMode != 0) 00966 { 00967 SG_INFO( " Time for setup: %.2lf\n", (float64_t)t/CLOCKS_PER_SEC); 00968 SG_INFO( 00969 "\n\n*********** Begin decomposition technique...\n"); 00970 } 00971 00972 /* arrays allocation */ 00973 bmem = SG_MALLOC(int32_t, ell); 00974 bmemrid = SG_MALLOC(int32_t, chunk_size); 00975 pbmr = SG_MALLOC(int32_t, chunk_size); 00976 cec = SG_MALLOC(int32_t, ell); 00977 indnzin = SG_MALLOC(int32_t, chunk_size); 00978 inold = SG_MALLOC(int32_t, chunk_size); 00979 incom = SG_MALLOC(int32_t, chunk_size); 00980 vau = SG_MALLOC(float64_t, ell); 00981 grad = SG_MALLOC(float64_t, ell); 00982 weight = SG_MALLOC(float64_t, dim); 00983 st = SG_MALLOC(float64_t, ell); 00984 stloc = SG_MALLOC(float64_t, ell); 00985 00986 for (i = 0; i < ell; i++) 00987 { 00988 bmem[i] = 0; 00989 cec[i] = 0; 00990 st[i] = 0; 00991 } 00992 00993 sp_y = SG_MALLOC(int32_t, chunk_size); 00994 sp_D = SG_MALLOC(float32_t, chunk_size*chunk_size); 00995 sp_alpha = SG_MALLOC(float64_t, chunk_size); 00996 sp_h = SG_MALLOC(float64_t, chunk_size); 00997 sp_hloc = SG_MALLOC(float64_t, chunk_size); 00998 00999 for (i = 0; i < chunk_size; i++) 01000 cec[index_in[i]] = cec[index_in[i]]+1; 01001 01002 for (i = chunk_size-1; i >= 0; i--) 01003 { 01004 incom[i] = -1; 01005 sp_alpha[i] = 0.0; 01006 bmem[index_in[i]] = 1; 01007 } 01008 01009 if (verbosity == 1) 01010 { 01011 SG_INFO( " IT | Prep Time | Solver IT | Solver Time |"); 01012 SG_INFO( " Grad Time | KKT violation\n"); 01013 SG_INFO( "------+-----------+-----------+-------------+"); 01014 SG_INFO( "-----------+--------------\n"); 01015 } 01016 01017 /***************************************************************************/ 01018 /* Begin the problem resolution loop */ 01019 nit = 0; 01020 do 01021 { 01022 t = clock(); 01023 if ((nit % 10) == 9) 01024 { 01025 if ((t-ti) > 0) 01026 total_time += (float64_t)(t-ti) / CLOCKS_PER_SEC; 01027 else 01028 total_time += (float64_t)(ti-t) / CLOCKS_PER_SEC; 01029 ti = t; 01030 } 01031 01032 if (verbosity > 1) 01033 SG_INFO("\n*********** ITERATION: %d\n", nit + 1); 01034 else if (verbosity > 0) 01035 SG_INFO( "%5d |", nit + 1); 01036 else 01037 SG_INFO( "."); 01038 fflush(stdout); 01039 01040 nzout = 0; 01041 for (k = 0; k < nb; k++) 01042 if (alpha_out(k)>DELTAsv) 01043 { 01044 indnzout[nzout] = index_out[k]; 01045 nzout++; 01046 } 01047 01048 sp_e = 0.; 01049 for (j = 0; j < nzout; j++) 01050 { 01051 vau[j] = y[indnzout[j]]*alpha[indnzout[j]]; 01052 sp_e += vau[j]; 01053 } 01054 01055 if (verbosity > 1) 01056 SG_INFO( " spe: %e ", sp_e); 01057 01058 for (i = 0; i < chunk_size; i++) 01059 sp_y[i] = y_in(i); 01060 01061 /* Construct the objective function Hessian */ 01062 for (i = 0; i < chunk_size; i++) 01063 { 01064 iin = index_in[i]; 01065 ptmw = Cache->GetRow(iin); 01066 if (ptmw != 0) 01067 { 01068 for (j = 0; j <= i; j++) 01069 sp_D[i*chunk_size + j] = sp_y[i]*sp_y[j] * ptmw[index_in[j]]; 01070 } 01071 else if (incom[i] == -1) 01072 for (j = 0; j <= i; j++) 01073 sp_D[i*chunk_size + j] = sp_y[i]*sp_y[j] 01074 * (float32_t)KER->Get(iin, index_in[j]); 01075 else 01076 { 01077 for (j = 0; j < i; j++) 01078 if (incom[j] == -1) 01079 sp_D[i*chunk_size + j] 01080 = sp_y[i]*sp_y[j] * (float32_t)KER->Get(iin, index_in[j]); 01081 else 01082 sp_D[i*chunk_size + j] 01083 = sp_D[incom[j]*chunk_size + incom[i]]; 01084 sp_D[i*chunk_size + i] 01085 = sp_y[i]*sp_y[i] * (float32_t)KER->Get(iin, index_in[i]); 01086 } 01087 } 01088 for (i = 0; i < chunk_size; i++) 01089 { 01090 for (j = 0; j < i; j++) 01091 sp_D[j*chunk_size + i] = sp_D[i*chunk_size + j]; 01092 } 01093 01094 if (nit == 0 && PreprocessMode > 0) 01095 { 01096 for (i = 0; i < chunk_size; i++) 01097 { 01098 iin = index_in[i]; 01099 aux = 0.; 01100 ptmw = Cache->GetRow(iin); 01101 if (ptmw == NULL) 01102 for (j = 0; j < nzout; j++) 01103 aux += vau[j] * KER->Get(iin, indnzout[j]); 01104 else 01105 for (j = 0; j < nzout; j++) 01106 aux += vau[j] * ptmw[indnzout[j]]; 01107 sp_h[i] = y_in(i) * aux - 1.0; 01108 } 01109 } 01110 else 01111 { 01112 for (i = 0; i < chunk_size; i++) 01113 vau[i] = alpha_in(i); 01114 for (i = 0; i < chunk_size; i++) 01115 { 01116 aux = 0.0; 01117 for (j = 0; j < chunk_size; j++) 01118 aux += sp_D[i*chunk_size + j] * vau[j]; 01119 sp_h[i] = st[index_in[i]] * y_in(i) - 1.0 - aux; 01120 } 01121 } 01122 01123 t = clock() - t; 01124 if (verbosity > 1) 01125 SG_INFO( 01126 " Preparation Time: %.2lf\n", (float64_t)t/CLOCKS_PER_SEC); 01127 else if (verbosity > 0) 01128 SG_INFO( " %8.2lf |", (float64_t)t/CLOCKS_PER_SEC); 01129 tot_prep_time += (float64_t)t/CLOCKS_PER_SEC; 01130 01131 /*** Proximal point modification: first type ***/ 01132 01133 if (tau_proximal < 0.0) 01134 tau_proximal_this = 0.0; 01135 else 01136 tau_proximal_this = tau_proximal; 01137 proximal_count = 0; 01138 do { 01139 t = clock(); 01140 for (i = 0; i < chunk_size; i++) 01141 { 01142 vau[i] = sp_D[i*chunk_size + i]; 01143 sp_h[i] -= tau_proximal_this * alpha_in(i); 01144 sp_D[i*chunk_size + i] += (float32_t)tau_proximal_this; 01145 } 01146 01147 if (kktold < delta*100) 01148 for (i = 0; i < chunk_size; i++) 01149 sp_alpha[i] = alpha_in(i); 01150 else 01151 for (i = 0; i < chunk_size; i++) 01152 sp_alpha[i] = 0.0; 01153 01154 /*** call the chosen inner gradient projection QP solver ***/ 01155 i = gpm_solver(projection_solver, projection_projector, chunk_size, 01156 sp_D, sp_h, c_const, sp_e, sp_y, sp_alpha, 01157 DELTAvpm*DELTAkin, &lsCount, &projCount); 01158 01159 if (i > vpmWarningThreshold) 01160 { 01161 if (ker_type == 2) 01162 { 01163 SG_INFO("\n WARNING: inner subproblem hard to solve;"); 01164 SG_INFO(" setting a smaller -q or"); 01165 SG_INFO(" tuning -c and -g options might help.\n"); 01166 } 01167 else 01168 { 01169 SG_INFO("\n WARNING: inner subproblem hard to solve;"); 01170 SG_INFO(" set a smaller -q or"); 01171 SG_INFO(" try a better data scaling.\n"); 01172 } 01173 } 01174 01175 t = clock() - t; 01176 tot_vpm_iter += i; 01177 tot_vpm_secant += projCount; 01178 tot_vpm_time += (float64_t)t/CLOCKS_PER_SEC; 01179 if (verbosity > 1) 01180 { 01181 SG_INFO(" Solver it: %d", i); 01182 SG_INFO(", ls: %d", lsCount); 01183 SG_INFO(", time: %.2lf\n", (float64_t)t/CLOCKS_PER_SEC); 01184 } 01185 else if (verbosity > 0) 01186 { 01187 SG_INFO(" %6d", i); 01188 SG_INFO(" | %8.2lf |", (float64_t)t/CLOCKS_PER_SEC); 01189 } 01190 01191 /*** Proximal point modification: second type ***/ 01192 01193 for (i = 0; i < chunk_size; i++) 01194 sp_D[i*chunk_size + i] = (float32_t)vau[i]; 01195 tau_proximal_this = 0.0; 01196 if (tau_proximal < 0.0) 01197 { 01198 dfval = 0.0; 01199 for (i = 0; i < chunk_size; i++) 01200 { 01201 aux = 0.0; 01202 for (j = 0; j < chunk_size; j++) 01203 aux += sp_D[i*chunk_size + j]*(alpha_in(j) - sp_alpha[j]); 01204 dfval += (0.5*aux - st[index_in[i]]*y_in(i) + 1.0) * (alpha_in(i) - sp_alpha[i]); 01205 } 01206 01207 aux=0.0; 01208 for (i = 0; i < chunk_size; i++) 01209 aux += (alpha_in(i) - sp_alpha[i])*(alpha_in(i) - sp_alpha[i]); 01210 01211 if ((-dfval/aux) < -0.5*tau_proximal) 01212 { 01213 tau_proximal_this = -tau_proximal; 01214 if (verbosity > 0) 01215 SG_DEBUG("tau threshold: %lf ", -dfval/aux); 01216 } 01217 } 01218 proximal_count++; 01219 } while (tau_proximal_this != 0.0 && proximal_count < 2); // Proximal point loop 01220 01221 t = clock(); 01222 01223 nzin = 0; 01224 for (j = 0; j < chunk_size; j++) 01225 { 01226 if (nit == 0) 01227 aux = sp_alpha[j]; 01228 else 01229 aux = sp_alpha[j] - alpha_in(j); 01230 if (fabs(aux) > DELTAsv) 01231 { 01232 indnzin[nzin] = index_in[j]; 01233 grad[nzin] = aux * y_in(j); 01234 nzin++; 01235 } 01236 } 01237 01238 // in case of LINADD enabled use faster linadd variant 01239 if (KER->get_kernel()->has_property(KP_LINADD) && get_linadd_enabled()) 01240 { 01241 KER->get_kernel()->clear_normal() ; 01242 01243 for (j = 0; j < nzin; j++) 01244 KER->get_kernel()->add_to_normal(indnzin[j], grad[j]); 01245 01246 if (nit == 0 && PreprocessMode > 0) 01247 { 01248 for (j = 0; j < nzout; j++) 01249 { 01250 jin = indnzout[j]; 01251 KER->get_kernel()->add_to_normal(jin, alpha[jin] * y[jin]); 01252 } 01253 } 01254 01255 for (i = 0; i < ell; i++) 01256 st[i] += KER->get_kernel()->compute_optimized(i); 01257 } 01258 else // nonlinear kernel 01259 { 01260 k = Cache->DivideMP(ing, indnzin, nzin); 01261 for (j = 0; j < k; j++) 01262 { 01263 ptmw = Cache->FillRow(indnzin[ing[j]]); 01264 for (i = 0; i < ell; i++) 01265 st[i] += grad[ing[j]] * ptmw[i]; 01266 } 01267 01268 if (PreprocessMode > 0 && nit == 0) 01269 { 01270 clock_t ti2; 01271 01272 ti2 = clock(); 01273 for (j = 0; j < nzout; j++) 01274 { 01275 jin = indnzout[j]; 01276 ptmw = Cache->FillRow(jin); 01277 for (i = 0; i < ell; i++) 01278 st[i] += alpha[jin] * y[jin] * ptmw[i]; 01279 } 01280 if (verbosity > 1) 01281 SG_INFO( 01282 " G*x0 time: %.2lf\n", (float64_t)(clock()-ti2)/CLOCKS_PER_SEC); 01283 } 01284 } 01285 01286 /*** sort the vectors for cache managing ***/ 01287 01288 t = clock() - t; 01289 if (verbosity > 1) 01290 SG_INFO( 01291 " Gradient updating time: %.2lf\n", (float64_t)t/CLOCKS_PER_SEC); 01292 else if (verbosity > 0) 01293 SG_INFO( " %8.2lf |", (float64_t)t/CLOCKS_PER_SEC); 01294 tot_st_time += (float64_t)t/CLOCKS_PER_SEC; 01295 01296 /*** global updating of the solution vector ***/ 01297 for (i = 0; i < chunk_size; i++) 01298 alpha_in(i) = sp_alpha[i]; 01299 01300 if (verbosity > 1) 01301 { 01302 j = k = 0; 01303 for (i = 0; i < ell; i++) 01304 { 01305 if (is_free(i)) j++; 01306 if (is_bound(i)) k++; 01307 } 01308 SG_INFO(" SV: %d", j+k); 01309 SG_INFO(", BSV: %d\n", k); 01310 } 01311 Cache->Iteration(); 01312 nit = nit+1; 01313 } while (!optimal() && !(CSignal::cancel_computations())); 01314 /* End of the problem resolution loop */ 01315 /***************************************************************************/ 01316 01317 t = clock(); 01318 if ((t-ti) > 0) 01319 total_time += (float64_t)(t-ti) / CLOCKS_PER_SEC; 01320 else 01321 total_time += (float64_t)(ti-t) / CLOCKS_PER_SEC; 01322 ti = t; 01323 01324 memcpy(solution, alpha, ell * sizeof(float64_t)); 01325 01326 /* objective function evaluation */ 01327 fval = 0.0; 01328 for (i = 0; i < ell; i++) 01329 fval += alpha[i]*(y[i]*st[i]*0.5 - 1.0); 01330 01331 SG_INFO("\n------+-----------+-----------+-------------+"); 01332 SG_INFO("-----------+--------------\n"); 01333 SG_INFO( 01334 "\n- TOTAL ITERATIONS: %i\n", nit); 01335 01336 if (verbosity > 1) 01337 { 01338 j = 0; 01339 k = 0; 01340 z = 0; 01341 for (i = ell - 1; i >= 0; i--) 01342 { 01343 if (cec[i] > 0) j++; 01344 if (cec[i] > 1) k++; 01345 if (cec[i] > 2) z++; 01346 } 01347 SG_INFO( 01348 "- Variables entering the working set at least one time: %i\n", j); 01349 SG_INFO( 01350 "- Variables entering the working set at least two times: %i\n", k); 01351 SG_INFO( 01352 "- Variables entering the working set at least three times: %i\n", z); 01353 } 01354 01355 01356 SG_FREE(bmem); 01357 SG_FREE(bmemrid); 01358 SG_FREE(pbmr); 01359 SG_FREE(cec); 01360 SG_FREE(ing); 01361 SG_FREE(inaux); 01362 SG_FREE(indnzin); 01363 SG_FREE(index_in); 01364 SG_FREE(inold); 01365 SG_FREE(incom); 01366 SG_FREE(indnzout); 01367 SG_FREE(index_out); 01368 SG_FREE(vau); 01369 SG_FREE(alpha); 01370 SG_FREE(weight); 01371 SG_FREE(grad); 01372 SG_FREE(stloc); 01373 SG_FREE(st); 01374 SG_FREE(sp_h); 01375 SG_FREE(sp_hloc); 01376 SG_FREE(sp_y); 01377 SG_FREE(sp_D); 01378 SG_FREE(sp_alpha); 01379 delete Cache; 01380 01381 aux = KER->KernelEvaluations; 01382 SG_INFO( "- Total CPU time: %lf\n", total_time); 01383 if (verbosity > 0) 01384 { 01385 SG_INFO( 01386 "- Total kernel evaluations: %.0lf\n", aux); 01387 SG_INFO( 01388 "- Total inner solver iterations: %i\n", tot_vpm_iter); 01389 if (projection_projector == 1) 01390 SG_INFO( 01391 "- Total projector iterations: %i\n", tot_vpm_secant); 01392 SG_INFO( 01393 "- Total preparation time: %lf\n", tot_prep_time); 01394 SG_INFO( 01395 "- Total inner solver time: %lf\n", tot_vpm_time); 01396 SG_INFO( 01397 "- Total gradient updating time: %lf\n", tot_st_time); 01398 } 01399 SG_INFO( "- Objective function value: %lf\n", fval); 01400 objective_value=fval; 01401 return aux; 01402 } 01403 01404 /******************************************************************************/ 01405 /*** Quick sort for integer vectors ***/ 01406 /******************************************************************************/ 01407 void quick_si(int32_t a[], int32_t n) 01408 { 01409 int32_t i, j, s, d, l, x, w, ps[20], pd[20]; 01410 01411 l = 0; 01412 ps[0] = 0; 01413 pd[0] = n-1; 01414 do 01415 { 01416 s = ps[l]; 01417 d = pd[l]; 01418 l--; 01419 do 01420 { 01421 i = s; 01422 j = d; 01423 x = a[(s+d)/2]; 01424 do 01425 { 01426 while (a[i] < x) i++; 01427 while (a[j] > x) j--; 01428 if (i <= j) 01429 { 01430 w = a[i]; 01431 a[i] = a[j]; 01432 i++; 01433 a[j] = w; 01434 j--; 01435 } 01436 } while(i<=j); 01437 if (j-s > d-i) 01438 { 01439 l++; 01440 ps[l] = s; 01441 pd[l] = j; 01442 s = i; 01443 } 01444 else 01445 { 01446 if (i < d) 01447 { 01448 l++; 01449 ps[l] = i; 01450 pd[l] = d; 01451 } 01452 d = j; 01453 } 01454 } while (s < d); 01455 } while (l >= 0); 01456 } 01457 01458 /******************************************************************************/ 01459 /*** Quick sort for real vectors returning also the exchanges ***/ 01460 /******************************************************************************/ 01461 void quick_s2(float64_t a[], int32_t n, int32_t ia[]) 01462 { 01463 int32_t i, j, s, d, l, iw, ps[20], pd[20]; 01464 float64_t x, w; 01465 01466 l = 0; 01467 ps[0] = 0; 01468 pd[0] = n-1; 01469 do 01470 { 01471 s = ps[l]; 01472 d = pd[l]; 01473 l--; 01474 do 01475 { 01476 i = s; 01477 j = d; 01478 x = a[(s+d)/2]; 01479 do 01480 { 01481 while (a[i] < x) i++; 01482 while (a[j] > x) j--; 01483 if (i <= j) 01484 { 01485 iw = ia[i]; 01486 w = a[i]; 01487 ia[i] = ia[j]; 01488 a[i] = a[j]; 01489 i++; 01490 ia[j] = iw; 01491 a[j] = w; 01492 j--; 01493 } 01494 } while (i <= j); 01495 if (j-s > d-i) 01496 { 01497 l++; 01498 ps[l] = s; 01499 pd[l] = j; 01500 s = i; 01501 } 01502 else 01503 { 01504 if (i < d) 01505 { 01506 l++; 01507 ps[l] = i; 01508 pd[l] = d; 01509 } 01510 d = j; 01511 } 01512 } while (s < d); 01513 } while(l>=0); 01514 } 01515 01516 /******************************************************************************/ 01517 /*** Quick sort for integer vectors returning also the exchanges ***/ 01518 /******************************************************************************/ 01519 void quick_s3(int32_t a[], int32_t n, int32_t ia[]) 01520 { 01521 int32_t i, j, s, d, l, iw, w, x, ps[20], pd[20]; 01522 01523 l = 0; 01524 ps[0] = 0; 01525 pd[0] = n-1; 01526 do 01527 { 01528 s = ps[l]; 01529 d = pd[l]; 01530 l--; 01531 do 01532 { 01533 i = s; 01534 j = d; 01535 x = a[(s+d)/2]; 01536 do 01537 { 01538 while (a[i] < x) i++; 01539 while (a[j] > x) j--; 01540 if (i <= j) 01541 { 01542 iw = ia[i]; 01543 w = a[i]; 01544 ia[i] = ia[j]; 01545 a[i] = a[j]; 01546 i++; 01547 ia[j] = iw; 01548 a[j] = w; 01549 j--; 01550 } 01551 } while (i <= j); 01552 if (j-s > d-i) 01553 { 01554 l++; 01555 ps[l] = s; 01556 pd[l] = j; 01557 s = i; 01558 } 01559 else 01560 { 01561 if (i < d) 01562 { 01563 l++; 01564 ps[l] = i; 01565 pd[l] = d; 01566 } 01567 d = j; 01568 } 01569 } while (s < d); 01570 } while (l >= 0); 01571 } 01572 } 01573 01574 #endif // DOXYGEN_SHOULD_SKIP_THIS 01575 01576 /******************************************************************************/ 01577 /*** End of gpdtsolve.cpp file ***/ 01578 /******************************************************************************/