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