00001
00002
00003
00004
00005
00006
00007
00008
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
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
00059 if (nev>n)
00060 SG_SERROR("Number of required eigenpairs is greater than order of the matrix\n");
00061
00062
00063 if (mode!=1 && mode!=3)
00064 SG_SERROR("Mode not supported yet\n");
00065
00066
00067
00068 int ido = 0;
00069
00070
00071
00072 char bmat[2] = "I";
00073 if (rhs)
00074 bmat[0] = 'G';
00075
00076
00077 double tol = tolerance;
00078
00079
00080 double* resid = SG_MALLOC(double, n);
00081
00082 for (i=0; i<n; i++)
00083 resid[i] = 1.0;
00084
00085
00086
00087 int ncv = nev*4>n ? n : nev*4;
00088
00089
00090 int ldv = n;
00091 double* v = SG_MALLOC(double, ldv*ncv);
00092
00093
00094 int* iparam = SG_MALLOC(int, 11);
00095
00096 iparam[0] = 1;
00097
00098 iparam[2] = 3*n;
00099
00100 iparam[6] = mode;
00101
00102
00103 int* ipntr = SG_CALLOC(int, 11);
00104
00105
00106 double* workd = SG_MALLOC(double, 3*n);
00107 int lworkl = ncv*(ncv+8);
00108 double* workl = SG_MALLOC(double, lworkl);
00109
00110
00111 int info = 1;
00112
00113
00114 char* which_ = strdup(which);
00115
00116 char* all_ = strdup("A");
00117
00118
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
00149 if (mode==3)
00150 {
00151
00152 if (shift!=0.0)
00153 {
00154 SG_SDEBUG("ARPACK: Subtracting shift\n");
00155
00156
00157 if (rhs && is_rhs_diag)
00158
00159 for (i=0; i<n; i++)
00160 matrix[i*n+i] -= rhs[i]*shift;
00161
00162 else
00163
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
00174 for (i=0; i<n*n; i++)
00175 {
00176 if (matrix[i]!=0.0)
00177 nnz++;
00178 }
00179
00180 double* val = doubleMalloc(nnz);
00181 int* rowind = intMalloc(nnz);
00182 int* colptr = intMalloc(n+1);
00183 nnz = 0;
00184
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
00200 dCreate_CompCol_Matrix(&slu_A,n,n,nnz,val,rowind,colptr,SLU_NC,SLU_D,SLU_GE);
00201
00202
00203 set_default_options(&options);
00204 options.Equil = YES;
00205 options.SymmetricMode = YES;
00206 StatInit(&stat);
00207
00208
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
00230 if (pos)
00231 {
00232
00233 SG_SDEBUG("ARPACK: Using Cholesky factorization.\n");
00234 clapack_dpotrf(CblasColMajor,CblasUpper,n,matrix,n);
00235 }
00236 else
00237 {
00238
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
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
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
00271 for (i=0; i<n; i++)
00272 slu_Bv[i] = (workd+ipntr[0]-1)[i];
00273 slu_info = 0;
00274
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
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
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
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
00412 int* select = SG_MALLOC(int, ncv);
00413
00414 double* d = SG_MALLOC(double, 2*ncv);
00415
00416 double sigma = shift;
00417
00418
00419 int ierr = 0;
00420
00421
00422 int rvec = 1;
00423
00424 SG_SDEBUG("APRACK: Starting DSEUPD.\n");
00425
00426
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
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
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
00452 SG_FREE(select);
00453 SG_FREE(d);
00454 }
00455
00456
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
00468 #endif