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