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 * Written (W) 2011 Sergey Lisitsyn 00008 * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society 00009 */ 00010 00011 #include <shogun/mathematics/arpack.h> 00012 #ifdef HAVE_ARPACK 00013 #ifdef HAVE_LAPACK 00014 #include <shogun/io/SGIO.h> 00015 #include <string.h> 00016 #include <shogun/mathematics/lapack.h> 00017 00018 #ifdef HAVE_SUPERLU 00019 #include <superlu/slu_ddefs.h> 00020 #endif 00021 00022 00024 extern "C" void dsaupd_(int *ido, char *bmat, int *n, char *which, 00025 int *nev, double *tol, double *resid, int *ncv, 00026 double *v, int *ldv, int *iparam, int *ipntr, 00027 double *workd, double *workl, int *lworkl, 00028 int *info); 00029 00031 extern "C" void dseupd_(int *rvec, char *All, int *select, double *d, 00032 double *v, int *ldv, double *sigma, 00033 char *bmat, int *n, char *which, int *nev, 00034 double *tol, double *resid, int *ncv, double *tv, 00035 int *tldv, int *iparam, int *ipntr, double *workd, 00036 double *workl, int *lworkl, int *ierr); 00037 00038 using namespace shogun; 00039 00040 namespace shogun 00041 { 00042 void arpack_dsxupd(double* matrix, double* rhs, bool is_rhs_diag, int n, int nev, 00043 const char* which, bool use_superlu, int mode, bool pos, 00044 double shift, double tolerance, double* eigenvalues, 00045 double* eigenvectors, int& status) 00046 { 00047 // loop vars 00048 int i,j; 00049 00050 if (use_superlu) 00051 { 00052 #ifndef HAVE_SUPERLU 00053 use_superlu=false; 00054 SG_SINFO("Required SUPERLU isn't available in this configuration\n"); 00055 #endif 00056 } 00057 00058 // check if nev is greater than n 00059 if (nev>n) 00060 SG_SERROR("Number of required eigenpairs is greater than order of the matrix\n"); 00061 00062 // check specified mode 00063 if (mode!=1 && mode!=3) 00064 SG_SERROR("Mode not supported yet\n"); 00065 00066 // init ARPACK's reverse communication parameter 00067 // (should be zero initially) 00068 int ido = 0; 00069 00070 // specify general or non-general eigenproblem will be solved 00071 // w.r.t. to given rhs_diag 00072 char bmat[2] = "I"; 00073 if (rhs) 00074 bmat[0] = 'G'; 00075 00076 // init tolerance (zero means machine precision) 00077 double tol = tolerance; 00078 00079 // allocate array to hold residuals 00080 double* resid = SG_MALLOC(double, n); 00081 // fill residual vector with ones 00082 for (i=0; i<n; i++) 00083 resid[i] = 1.0; 00084 00085 // set number of Lanczos basis vectors to be used 00086 // (with max(4*nev,n) sufficient for most tasks) 00087 int ncv = nev*4>n ? n : nev*4; 00088 00089 // allocate array 'v' for dsaupd routine usage 00090 int ldv = n; 00091 double* v = SG_MALLOC(double, ldv*ncv); 00092 00093 // init array for i/o params for routine 00094 int* iparam = SG_MALLOC(int, 11); 00095 // specify method for selecting implicit shifts (1 - exact shifts) 00096 iparam[0] = 1; 00097 // specify max number of iterations 00098 iparam[2] = 3*n; 00099 // set the computation mode (1 for regular or 3 for shift-inverse) 00100 iparam[6] = mode; 00101 00102 // init array indicating locations of vectors for routine callback 00103 int* ipntr = SG_CALLOC(int, 11); 00104 00105 // allocate workaround arrays 00106 double* workd = SG_MALLOC(double, 3*n); 00107 int lworkl = ncv*(ncv+8); 00108 double* workl = SG_MALLOC(double, lworkl); 00109 00110 // init info holding status (1 means that residual vector is provided initially) 00111 int info = 1; 00112 00113 // which eigenpairs to find 00114 char* which_ = strdup(which); 00115 // All 00116 char* all_ = strdup("A"); 00117 00118 // ipiv for LUP factorization 00119 int* ipiv = NULL; 00120 00121 #ifdef HAVE_SUPERLU 00122 char equed[1]; 00123 void* work = NULL; 00124 int lwork = 0; 00125 SuperMatrix slu_A, slu_L, slu_U, slu_B, slu_X; 00126 superlu_options_t options; 00127 SuperLUStat_t stat; 00128 mem_usage_t mem_usage; 00129 int *perm_c=NULL, *perm_r=NULL, *etree=NULL; 00130 double *R=NULL, *C=NULL; 00131 if (mode==3 && use_superlu) 00132 { 00133 perm_c = intMalloc(n); 00134 perm_r = intMalloc(n); 00135 etree = intMalloc(n); 00136 R = doubleMalloc(n); 00137 C = doubleMalloc(n); 00138 } 00139 double ferr; 00140 double berr; 00141 double rcond; 00142 double rpg; 00143 int slu_info; 00144 double* slu_Bv=NULL; 00145 double* slu_Xv=NULL; 00146 #endif 00147 00148 // shift-invert mode init 00149 if (mode==3) 00150 { 00151 // subtract shift from main diagonal if necessary 00152 if (shift!=0.0) 00153 { 00154 SG_SDEBUG("ARPACK: Subtracting shift\n"); 00155 // if right hand side diagonal matrix is provided 00156 00157 if (rhs && is_rhs_diag) 00158 // subtract I*diag(rhs_diag) 00159 for (i=0; i<n; i++) 00160 matrix[i*n+i] -= rhs[i]*shift; 00161 00162 else 00163 // subtract I 00164 for (i=0; i<n; i++) 00165 matrix[i*n+i] -= shift; 00166 } 00167 00168 if (use_superlu) 00169 { 00170 #ifdef HAVE_SUPERLU 00171 SG_SDEBUG("SUPERLU: Constructing sparse matrix.\n"); 00172 int nnz = 0; 00173 // get number of non-zero elements 00174 for (i=0; i<n*n; i++) 00175 { 00176 if (matrix[i]!=0.0) 00177 nnz++; 00178 } 00179 // allocate arrays to store sparse matrix 00180 double* val = doubleMalloc(nnz); 00181 int* rowind = intMalloc(nnz); 00182 int* colptr = intMalloc(n+1); 00183 nnz = 0; 00184 // construct sparse matrix 00185 for (i=0; i<n; i++) 00186 { 00187 colptr[i] = nnz; 00188 for (j=0; j<n; j++) 00189 { 00190 if (matrix[i*n+j]!=0.0) 00191 { 00192 val[nnz] = matrix[i*n+j]; 00193 rowind[nnz] = j; 00194 nnz++; 00195 } 00196 } 00197 } 00198 colptr[i] = nnz; 00199 // create CCS matrix 00200 dCreate_CompCol_Matrix(&slu_A,n,n,nnz,val,rowind,colptr,SLU_NC,SLU_D,SLU_GE); 00201 00202 // initialize options 00203 set_default_options(&options); 00204 options.Equil = YES; 00205 options.SymmetricMode = YES; 00206 StatInit(&stat); 00207 00208 // factorize 00209 slu_info = 0; 00210 slu_Bv = doubleMalloc(n); 00211 slu_Xv = doubleMalloc(n); 00212 dCreate_Dense_Matrix(&slu_B,n,1,slu_Bv,n,SLU_DN,SLU_D,SLU_GE); 00213 dCreate_Dense_Matrix(&slu_X,n,1,slu_Xv,n,SLU_DN,SLU_D,SLU_GE); 00214 slu_B.ncol = 0; 00215 SG_SDEBUG("SUPERLU: Factorizing\n"); 00216 dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C, 00217 &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond, &ferr, &berr, 00218 &mem_usage,&stat,&slu_info); 00219 slu_B.ncol = 1; 00220 if (slu_info) 00221 { 00222 SG_SERROR("SUPERLU: Failed to factorize matrix, got %d code\n", slu_info); 00223 } 00224 options.Fact = FACTORED; 00225 #endif 00226 } 00227 else 00228 { 00229 // compute factorization according to pos value 00230 if (pos) 00231 { 00232 // with Cholesky 00233 SG_SDEBUG("ARPACK: Using Cholesky factorization.\n"); 00234 clapack_dpotrf(CblasColMajor,CblasUpper,n,matrix,n); 00235 } 00236 else 00237 { 00238 // with LUP 00239 SG_SDEBUG("ARPACK: Using LUP factorization.\n"); 00240 ipiv = SG_CALLOC(int, n); 00241 clapack_dgetrf(CblasColMajor,n,n,matrix,n,ipiv); 00242 } 00243 } 00244 } 00245 // main computation loop 00246 SG_SDEBUG("ARPACK: Starting main computation DSAUPD loop.\n"); 00247 do 00248 { 00249 dsaupd_(&ido, bmat, &n, which_, &nev, &tol, resid, 00250 &ncv, v, &ldv, iparam, ipntr, workd, workl, 00251 &lworkl, &info); 00252 00253 if ((ido==1)||(ido==-1)||(ido==2)) 00254 { 00255 if (mode==1) 00256 { 00257 // compute (workd+ipntr[1]-1) = A*(workd+ipntr[0]-1) 00258 cblas_dsymv(CblasColMajor,CblasUpper, 00259 n,1.0,matrix,n, 00260 (workd+ipntr[0]-1),1, 00261 0.0,(workd+ipntr[1]-1),1); 00262 } 00263 if (mode==3) 00264 { 00265 if (!rhs) 00266 { 00267 if (use_superlu) 00268 { 00269 #ifdef HAVE_SUPERLU 00270 // treat workd+ipntr(0) as B 00271 for (i=0; i<n; i++) 00272 slu_Bv[i] = (workd+ipntr[0]-1)[i]; 00273 slu_info = 0; 00274 // solve 00275 dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C, 00276 &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond, 00277 &ferr, &berr, &mem_usage, &stat, &slu_info); 00278 if (slu_info) 00279 SG_SERROR("SUPERLU: GOT %d\n", slu_info); 00280 // move elements from resulting X to workd+ipntr(1) 00281 for (i=0; i<n; i++) 00282 (workd+ipntr[1]-1)[i] = slu_Xv[i]; 00283 #endif 00284 } 00285 else 00286 { 00287 for (i=0; i<n; i++) 00288 (workd+ipntr[1]-1)[i] = (workd+ipntr[0]-1)[i]; 00289 if (pos) 00290 clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n); 00291 else 00292 clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n); 00293 } 00294 } 00295 else 00296 // HAVE RHS 00297 { 00298 if (ido==-1) 00299 { 00300 if (use_superlu) 00301 { 00302 #ifdef HAVE_SUPERLU 00303 for (i=0; i<n; i++) 00304 slu_Bv[i] = (workd+ipntr[0]-1)[i]; 00305 slu_info = 0; 00306 dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C, 00307 &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond, 00308 &ferr, &berr, &mem_usage, &stat, &slu_info); 00309 if (slu_info) 00310 SG_SERROR("SUPERLU: GOT %d\n", slu_info); 00311 for (i=0; i<n; i++) 00312 (workd+ipntr[1]-1)[i] = slu_Xv[i]; 00313 #endif 00314 } 00315 else 00316 { 00317 for (i=0; i<n; i++) 00318 (workd+ipntr[1]-1)[i] = rhs[i]*(workd+ipntr[0]-1)[i]; 00319 if (pos) 00320 clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n); 00321 else 00322 clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n); 00323 } 00324 } 00325 if (ido==1) 00326 { 00327 if (use_superlu) 00328 { 00329 #ifdef HAVE_SUPERLU 00330 for (i=0; i<n; i++) 00331 slu_Bv[i] = (workd+ipntr[2]-1)[i]; 00332 slu_info = 0; 00333 dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C, 00334 &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond, 00335 &ferr, &berr, &mem_usage, &stat, &slu_info); 00336 if (slu_info) 00337 SG_SERROR("SUPERLU: GOT %d\n", slu_info); 00338 for (i=0; i<n; i++) 00339 (workd+ipntr[1]-1)[i] = slu_Xv[i]; 00340 #endif 00341 } 00342 else 00343 { 00344 for (i=0; i<n; i++) 00345 (workd+ipntr[1]-1)[i] = (workd+ipntr[2]-1)[i]; 00346 if (pos) 00347 clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n); 00348 else 00349 clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n); 00350 } 00351 } 00352 if (ido==2) 00353 { 00354 if (is_rhs_diag) 00355 { 00356 for (i=0; i<n; i++) 00357 (workd+ipntr[1]-1)[i] = rhs[i]*(workd+ipntr[0]-1)[i]; 00358 } 00359 else 00360 { 00361 cblas_dsymv(CblasColMajor,CblasUpper, 00362 n,1.0,rhs,n, 00363 (workd+ipntr[0]-1),1, 00364 0.0,(workd+ipntr[1]-1),1); 00365 } 00366 } 00367 } 00368 } 00369 } 00370 } while ((ido==1)||(ido==-1)||(ido==2)); 00371 00372 if (!pos && mode==3) SG_FREE(ipiv); 00373 00374 if (mode==3 && use_superlu) 00375 { 00376 #ifdef HAVE_SUPERLU 00377 SUPERLU_FREE(slu_Bv); 00378 SUPERLU_FREE(slu_Xv); 00379 SUPERLU_FREE(perm_r); 00380 SUPERLU_FREE(perm_c); 00381 SUPERLU_FREE(R); 00382 SUPERLU_FREE(C); 00383 SUPERLU_FREE(etree); 00384 if (lwork!=0) SUPERLU_FREE(work); 00385 Destroy_CompCol_Matrix(&slu_A); 00386 StatFree(&stat); 00387 Destroy_SuperMatrix_Store(&slu_B); 00388 Destroy_SuperMatrix_Store(&slu_X); 00389 Destroy_SuperNode_Matrix(&slu_L); 00390 Destroy_CompCol_Matrix(&slu_U); 00391 #endif 00392 } 00393 00394 // check if DSAUPD failed 00395 if (info<0) 00396 { 00397 if ((info<=-1)&&(info>=-6)) 00398 SG_SWARNING("ARPACK: DSAUPD failed. Wrong parameter passed.\n"); 00399 else if (info==-7) 00400 SG_SWARNING("ARPACK: DSAUPD failed. Workaround array size is not sufficient.\n"); 00401 else 00402 SG_SWARNING("ARPACK: DSAUPD failed. Error code: %d.\n", info); 00403 00404 status = info; 00405 } 00406 else 00407 { 00408 if (info==1) 00409 SG_SWARNING("ARPACK: Maximum number of iterations reached.\n"); 00410 00411 // allocate select for dseupd 00412 int* select = SG_MALLOC(int, ncv); 00413 // allocate d to hold eigenvalues 00414 double* d = SG_MALLOC(double, 2*ncv); 00415 // sigma for dseupd 00416 double sigma = shift; 00417 00418 // init ierr indicating dseupd possible errors 00419 int ierr = 0; 00420 00421 // specify that eigenvectors are going to be computed too 00422 int rvec = 1; 00423 00424 SG_SDEBUG("APRACK: Starting DSEUPD.\n"); 00425 00426 // call dseupd_ routine 00427 dseupd_(&rvec, all_, select, d, v, &ldv, &sigma, bmat, 00428 &n, which_, &nev, &tol, resid, &ncv, v, &ldv, 00429 iparam, ipntr, workd, workl, &lworkl, &ierr); 00430 00431 // check for errors 00432 if (ierr!=0) 00433 { 00434 SG_SWARNING("ARPACK: DSEUPD failed with status %d.\n", ierr); 00435 status = -1; 00436 } 00437 else 00438 { 00439 SG_SDEBUG("ARPACK: Storing eigenpairs.\n"); 00440 00441 // store eigenpairs to specified arrays 00442 for (i=0; i<nev; i++) 00443 { 00444 eigenvalues[i] = d[i]; 00445 00446 for (j=0; j<n; j++) 00447 eigenvectors[j*nev+i] = v[i*n+j]; 00448 } 00449 } 00450 00451 // cleanup 00452 SG_FREE(select); 00453 SG_FREE(d); 00454 } 00455 00456 // cleanup 00457 SG_FREE(all_); 00458 SG_FREE(which_); 00459 SG_FREE(resid); 00460 SG_FREE(v); 00461 SG_FREE(iparam); 00462 SG_FREE(ipntr); 00463 SG_FREE(workd); 00464 SG_FREE(workl); 00465 }; 00466 } 00467 #endif /* HAVE_LAPACK */ 00468 #endif /* HAVE_ARPACK */