19 #include <superlu/slu_ddefs.h>
24 extern "C" void dsaupd_(
int *ido,
char *bmat,
int *n,
char *which,
25 int *nev,
double *tol,
double *resid,
int *ncv,
26 double *v,
int *ldv,
int *iparam,
int *ipntr,
27 double *workd,
double *workl,
int *lworkl,
31 extern "C" void dseupd_(
int *rvec,
char *All,
int *select,
double *d,
32 double *v,
int *ldv,
double *sigma,
33 char *bmat,
int *n,
char *which,
int *nev,
34 double *tol,
double *resid,
int *ncv,
double *tv,
35 int *tldv,
int *iparam,
int *ipntr,
double *workd,
36 double *workl,
int *lworkl,
int *ierr);
38 using namespace shogun;
42 void arpack_dsxupd(
double* matrix,
double* rhs,
bool is_rhs_diag,
int n,
int nev,
43 const char* which,
bool use_superlu,
int mode,
bool pos,
bool cov,
44 double shift,
double tolerance,
double* eigenvalues,
45 double* eigenvectors,
int& status)
54 SG_SINFO(
"Required SUPERLU isn't available in this configuration\n");
60 SG_SERROR(
"Number of required eigenpairs is greater than order of the matrix\n");
63 if (mode!=1 && mode!=3)
77 double tol = tolerance;
87 int ncv = nev*4>n ? n : nev*4;
107 int lworkl = ncv*(ncv+8);
108 double* workl =
SG_MALLOC(
double, lworkl);
114 char* which_ = strdup(which);
116 char* all_ = strdup(
"A");
125 SuperMatrix slu_A, slu_L, slu_U, slu_B, slu_X;
126 superlu_options_t options;
128 mem_usage_t mem_usage;
129 int *perm_c=NULL, *perm_r=NULL, *etree=NULL;
130 double *R=NULL, *C=NULL;
131 if (mode==3 && use_superlu)
133 perm_c = intMalloc(n);
134 perm_r = intMalloc(n);
135 etree = intMalloc(n);
154 SG_SDEBUG(
"ARPACK: Subtracting shift\n");
157 if (rhs && is_rhs_diag)
160 matrix[i*n+i] -= rhs[i]*shift;
165 matrix[i*n+i] -= shift;
171 SG_SDEBUG(
"SUPERLU: Constructing sparse matrix.\n");
174 for (i=0; i<n*n; i++)
180 double* val = doubleMalloc(nnz);
181 int* rowind = intMalloc(nnz);
182 int* colptr = intMalloc(n+1);
190 if (matrix[i*n+j]!=0.0)
192 val[nnz] = matrix[i*n+j];
200 dCreate_CompCol_Matrix(&slu_A,n,n,nnz,val,rowind,colptr,SLU_NC,SLU_D,SLU_GE);
203 set_default_options(&options);
207 options.SymmetricMode = YES;
208 options.ColPerm = MMD_AT_PLUS_A;
209 options.DiagPivotThresh = 0.001;
214 slu_Bv = doubleMalloc(n);
215 slu_Xv = doubleMalloc(n);
216 dCreate_Dense_Matrix(&slu_B,n,1,slu_Bv,n,SLU_DN,SLU_D,SLU_GE);
217 dCreate_Dense_Matrix(&slu_X,n,1,slu_Xv,n,SLU_DN,SLU_D,SLU_GE);
220 dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
221 &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond, &ferr, &berr,
222 &mem_usage,&stat,&slu_info);
226 SG_SERROR(
"SUPERLU: Failed to factorize matrix, got %d code\n", slu_info);
228 options.Fact = FACTORED;
237 SG_SDEBUG(
"ARPACK: Using Cholesky factorization.\n");
238 clapack_dpotrf(CblasColMajor,CblasUpper,n,matrix,n);
243 SG_SDEBUG(
"ARPACK: Using LUP factorization.\n");
245 clapack_dgetrf(CblasColMajor,n,n,matrix,n,ipiv);
250 SG_SDEBUG(
"ARPACK: Starting main computation DSAUPD loop.\n");
253 dsaupd_(&ido, bmat, &n, which_, &nev, &tol, resid,
254 &ncv, v, &ldv, iparam, ipntr, workd, workl,
257 if ((ido==1)||(ido==-1)||(ido==2))
264 cblas_dsymv(CblasColMajor,CblasUpper,
266 (workd+ipntr[0]-1),1,
267 0.0,(workd+ipntr[1]-1),1);
271 cblas_dsymv(CblasColMajor,CblasUpper,
273 (workd+ipntr[0]-1),1,
274 0.0,(workd+ipntr[1]-1),1);
275 cblas_dsymv(CblasColMajor,CblasUpper,
277 (workd+ipntr[1]-1),1,
278 0.0,(workd+ipntr[0]-1),1);
279 cblas_dcopy(n,workd+ipntr[0]-1,1,workd+ipntr[1]-1,1);
291 slu_Bv[i] = (workd+ipntr[0]-1)[i];
294 dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
295 &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond,
296 &ferr, &berr, &mem_usage, &stat, &slu_info);
298 SG_SERROR(
"SUPERLU: GOT %d\n", slu_info);
301 (workd+ipntr[1]-1)[i] = slu_Xv[i];
307 (workd+ipntr[1]-1)[i] = (workd+ipntr[0]-1)[i];
309 clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
311 clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
323 slu_Bv[i] = (workd+ipntr[0]-1)[i];
325 dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
326 &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond,
327 &ferr, &berr, &mem_usage, &stat, &slu_info);
329 SG_SERROR(
"SUPERLU: GOT %d\n", slu_info);
331 (workd+ipntr[1]-1)[i] = slu_Xv[i];
337 (workd+ipntr[1]-1)[i] = rhs[i]*(workd+ipntr[0]-1)[i];
339 clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
341 clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
350 slu_Bv[i] = (workd+ipntr[2]-1)[i];
352 dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
353 &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond,
354 &ferr, &berr, &mem_usage, &stat, &slu_info);
356 SG_SERROR(
"SUPERLU: GOT %d\n", slu_info);
358 (workd+ipntr[1]-1)[i] = slu_Xv[i];
364 (workd+ipntr[1]-1)[i] = (workd+ipntr[2]-1)[i];
366 clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
368 clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
376 (workd+ipntr[1]-1)[i] = rhs[i]*(workd+ipntr[0]-1)[i];
380 cblas_dsymv(CblasColMajor,CblasUpper,
382 (workd+ipntr[0]-1),1,
383 0.0,(workd+ipntr[1]-1),1);
389 }
while ((ido==1)||(ido==-1)||(ido==2));
391 if (!pos && mode==3)
SG_FREE(ipiv);
393 if (mode==3 && use_superlu)
396 SUPERLU_FREE(slu_Bv);
397 SUPERLU_FREE(slu_Xv);
398 SUPERLU_FREE(perm_r);
399 SUPERLU_FREE(perm_c);
403 if (lwork!=0) SUPERLU_FREE(work);
404 Destroy_CompCol_Matrix(&slu_A);
406 Destroy_SuperMatrix_Store(&slu_B);
407 Destroy_SuperMatrix_Store(&slu_X);
408 Destroy_SuperNode_Matrix(&slu_L);
409 Destroy_CompCol_Matrix(&slu_U);
416 if ((info<=-1)&&(info>=-6))
417 SG_SWARNING(
"ARPACK: DSAUPD failed. Wrong parameter passed.\n");
419 SG_SWARNING(
"ARPACK: DSAUPD failed. Workaround array size is not sufficient.\n");
421 SG_SWARNING(
"ARPACK: DSAUPD failed. Error code: %d.\n", info);
428 SG_SWARNING(
"ARPACK: Maximum number of iterations reached.\n");
435 double sigma = shift;
446 dseupd_(&rvec, all_, select, d, v, &ldv, &sigma, bmat,
447 &n, which_, &nev, &tol, resid, &ncv, v, &ldv,
448 iparam, ipntr, workd, workl, &lworkl, &ierr);
453 SG_SWARNING(
"ARPACK: DSEUPD failed with status %d.\n", ierr);
458 SG_SDEBUG(
"ARPACK: Storing eigenpairs.\n");
461 for (i=0; i<nev; i++)
463 eigenvalues[i] = d[i];
466 eigenvectors[j*nev+i] = v[i*n+j];