SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
arpack.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2011 Sergey Lisitsyn
8  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
9  */
10 
12 #ifdef HAVE_ARPACK
13 #ifdef HAVE_LAPACK
14 #include <shogun/io/SGIO.h>
15 #include <string.h>
17 
18 #ifdef HAVE_SUPERLU
19 #include <superlu/slu_ddefs.h>
20 #endif
21 
22 
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,
28  int *info);
29 
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);
37 
38 using namespace shogun;
39 
40 namespace shogun
41 {
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)
46 {
47  // loop vars
48  int i,j;
49 
50  if (use_superlu)
51  {
52 #ifndef HAVE_SUPERLU
53  use_superlu=false;
54  SG_SINFO("Required SUPERLU isn't available in this configuration\n");
55 #endif
56  }
57 
58  // check if nev is greater than n
59  if (nev>n)
60  SG_SERROR("Number of required eigenpairs is greater than order of the matrix\n");
61 
62  // check specified mode
63  if (mode!=1 && mode!=3)
64  SG_SERROR("Mode not supported yet\n");
65 
66  // init ARPACK's reverse communication parameter
67  // (should be zero initially)
68  int ido = 0;
69 
70  // specify general or non-general eigenproblem will be solved
71  // w.r.t. to given rhs_diag
72  char bmat[2] = "I";
73  if (rhs)
74  bmat[0] = 'G';
75 
76  // init tolerance (zero means machine precision)
77  double tol = tolerance;
78 
79  // allocate array to hold residuals
80  double* resid = SG_MALLOC(double, n);
81  // fill residual vector with ones
82  for (i=0; i<n; i++)
83  resid[i] = 1.0;
84 
85  // set number of Lanczos basis vectors to be used
86  // (with max(4*nev,n) sufficient for most tasks)
87  int ncv = nev*4>n ? n : nev*4;
88 
89  // allocate array 'v' for dsaupd routine usage
90  int ldv = n;
91  double* v = SG_MALLOC(double, ldv*ncv);
92 
93  // init array for i/o params for routine
94  int* iparam = SG_MALLOC(int, 11);
95  // specify method for selecting implicit shifts (1 - exact shifts)
96  iparam[0] = 1;
97  // specify max number of iterations
98  iparam[2] = 3*n;
99  // set the computation mode (1 for regular or 3 for shift-inverse)
100  iparam[6] = mode;
101 
102  // init array indicating locations of vectors for routine callback
103  int* ipntr = SG_CALLOC(int, 11);
104 
105  // allocate workaround arrays
106  double* workd = SG_MALLOC(double, 3*n);
107  int lworkl = ncv*(ncv+8);
108  double* workl = SG_MALLOC(double, lworkl);
109 
110  // init info holding status (1 means that residual vector is provided initially)
111  int info = 0;
112 
113  // which eigenpairs to find
114  char* which_ = strdup(which);
115  // All
116  char* all_ = strdup("A");
117 
118  // ipiv for LUP factorization
119  int* ipiv = NULL;
120 
121 #ifdef HAVE_SUPERLU
122  char equed[1];
123  void* work = NULL;
124  int lwork = 0;
125  SuperMatrix slu_A, slu_L, slu_U, slu_B, slu_X;
126  superlu_options_t options;
127  SuperLUStat_t stat;
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)
132  {
133  perm_c = intMalloc(n);
134  perm_r = intMalloc(n);
135  etree = intMalloc(n);
136  R = doubleMalloc(n);
137  C = doubleMalloc(n);
138  }
139  double ferr;
140  double berr;
141  double rcond;
142  double rpg;
143  int slu_info;
144  double* slu_Bv=NULL;
145  double* slu_Xv=NULL;
146 #endif
147 
148  // shift-invert mode init
149  if (mode==3)
150  {
151  // subtract shift from main diagonal if necessary
152  if (shift!=0.0)
153  {
154  SG_SDEBUG("ARPACK: Subtracting shift\n");
155  // if right hand side diagonal matrix is provided
156 
157  if (rhs && is_rhs_diag)
158  // subtract I*diag(rhs_diag)
159  for (i=0; i<n; i++)
160  matrix[i*n+i] -= rhs[i]*shift;
161 
162  else
163  // subtract I
164  for (i=0; i<n; i++)
165  matrix[i*n+i] -= shift;
166  }
167 
168  if (use_superlu)
169  {
170 #ifdef HAVE_SUPERLU
171  SG_SDEBUG("SUPERLU: Constructing sparse matrix.\n");
172  int nnz = 0;
173  // get number of non-zero elements
174  for (i=0; i<n*n; i++)
175  {
176  if (matrix[i]!=0.0)
177  nnz++;
178  }
179  // allocate arrays to store sparse matrix
180  double* val = doubleMalloc(nnz);
181  int* rowind = intMalloc(nnz);
182  int* colptr = intMalloc(n+1);
183  nnz = 0;
184  // construct sparse matrix
185  for (i=0; i<n; i++)
186  {
187  colptr[i] = nnz;
188  for (j=0; j<n; j++)
189  {
190  if (matrix[i*n+j]!=0.0)
191  {
192  val[nnz] = matrix[i*n+j];
193  rowind[nnz] = j;
194  nnz++;
195  }
196  }
197  }
198  colptr[i] = nnz;
199  // create CCS matrix
200  dCreate_CompCol_Matrix(&slu_A,n,n,nnz,val,rowind,colptr,SLU_NC,SLU_D,SLU_GE);
201 
202  // initialize options
203  set_default_options(&options);
204  //options.Equil = YES;
205  //options.SymmetricMode = YES;
206 
207  options.SymmetricMode = YES;
208  options.ColPerm = MMD_AT_PLUS_A;
209  options.DiagPivotThresh = 0.001;
210  StatInit(&stat);
211 
212  // factorize
213  slu_info = 0;
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);
218  slu_B.ncol = 0;
219  SG_SDEBUG("SUPERLU: Factorizing\n");
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);
223  slu_B.ncol = 1;
224  if (slu_info)
225  {
226  SG_SERROR("SUPERLU: Failed to factorize matrix, got %d code\n", slu_info);
227  }
228  options.Fact = FACTORED;
229 #endif
230  }
231  else
232  {
233  // compute factorization according to pos value
234  if (pos)
235  {
236  // with Cholesky
237  SG_SDEBUG("ARPACK: Using Cholesky factorization.\n");
238  clapack_dpotrf(CblasColMajor,CblasUpper,n,matrix,n);
239  }
240  else
241  {
242  // with LUP
243  SG_SDEBUG("ARPACK: Using LUP factorization.\n");
244  ipiv = SG_CALLOC(int, n);
245  clapack_dgetrf(CblasColMajor,n,n,matrix,n,ipiv);
246  }
247  }
248  }
249  // main computation loop
250  SG_SDEBUG("ARPACK: Starting main computation DSAUPD loop.\n");
251  do
252  {
253  dsaupd_(&ido, bmat, &n, which_, &nev, &tol, resid,
254  &ncv, v, &ldv, iparam, ipntr, workd, workl,
255  &lworkl, &info);
256 
257  if ((ido==1)||(ido==-1)||(ido==2))
258  {
259  if (mode==1)
260  {
261  if (!cov)
262  {
263  // compute (workd+ipntr[1]-1) = A*(workd+ipntr[0]-1)
264  cblas_dsymv(CblasColMajor,CblasUpper,
265  n,1.0,matrix,n,
266  (workd+ipntr[0]-1),1,
267  0.0,(workd+ipntr[1]-1),1);
268  }
269  else
270  {
271  cblas_dsymv(CblasColMajor,CblasUpper,
272  n,1.0,matrix,n,
273  (workd+ipntr[0]-1),1,
274  0.0,(workd+ipntr[1]-1),1);
275  cblas_dsymv(CblasColMajor,CblasUpper,
276  n,1.0,matrix,n,
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);
280  }
281  }
282  if (mode==3)
283  {
284  if (!rhs)
285  {
286  if (use_superlu)
287  {
288 #ifdef HAVE_SUPERLU
289  // treat workd+ipntr(0) as B
290  for (i=0; i<n; i++)
291  slu_Bv[i] = (workd+ipntr[0]-1)[i];
292  slu_info = 0;
293  // solve
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);
297  if (slu_info)
298  SG_SERROR("SUPERLU: GOT %d\n", slu_info);
299  // move elements from resulting X to workd+ipntr(1)
300  for (i=0; i<n; i++)
301  (workd+ipntr[1]-1)[i] = slu_Xv[i];
302 #endif
303  }
304  else
305  {
306  for (i=0; i<n; i++)
307  (workd+ipntr[1]-1)[i] = (workd+ipntr[0]-1)[i];
308  if (pos)
309  clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
310  else
311  clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
312  }
313  }
314  else
315  // HAVE RHS
316  {
317  if (ido==-1)
318  {
319  if (use_superlu)
320  {
321 #ifdef HAVE_SUPERLU
322  for (i=0; i<n; i++)
323  slu_Bv[i] = (workd+ipntr[0]-1)[i];
324  slu_info = 0;
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);
328  if (slu_info)
329  SG_SERROR("SUPERLU: GOT %d\n", slu_info);
330  for (i=0; i<n; i++)
331  (workd+ipntr[1]-1)[i] = slu_Xv[i];
332 #endif
333  }
334  else
335  {
336  for (i=0; i<n; i++)
337  (workd+ipntr[1]-1)[i] = rhs[i]*(workd+ipntr[0]-1)[i];
338  if (pos)
339  clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
340  else
341  clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
342  }
343  }
344  if (ido==1)
345  {
346  if (use_superlu)
347  {
348 #ifdef HAVE_SUPERLU
349  for (i=0; i<n; i++)
350  slu_Bv[i] = (workd+ipntr[2]-1)[i];
351  slu_info = 0;
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);
355  if (slu_info)
356  SG_SERROR("SUPERLU: GOT %d\n", slu_info);
357  for (i=0; i<n; i++)
358  (workd+ipntr[1]-1)[i] = slu_Xv[i];
359 #endif
360  }
361  else
362  {
363  for (i=0; i<n; i++)
364  (workd+ipntr[1]-1)[i] = (workd+ipntr[2]-1)[i];
365  if (pos)
366  clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
367  else
368  clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
369  }
370  }
371  if (ido==2)
372  {
373  if (is_rhs_diag)
374  {
375  for (i=0; i<n; i++)
376  (workd+ipntr[1]-1)[i] = rhs[i]*(workd+ipntr[0]-1)[i];
377  }
378  else
379  {
380  cblas_dsymv(CblasColMajor,CblasUpper,
381  n,1.0,rhs,n,
382  (workd+ipntr[0]-1),1,
383  0.0,(workd+ipntr[1]-1),1);
384  }
385  }
386  }
387  }
388  }
389  } while ((ido==1)||(ido==-1)||(ido==2));
390 
391  if (!pos && mode==3) SG_FREE(ipiv);
392 
393  if (mode==3 && use_superlu)
394  {
395 #ifdef HAVE_SUPERLU
396  SUPERLU_FREE(slu_Bv);
397  SUPERLU_FREE(slu_Xv);
398  SUPERLU_FREE(perm_r);
399  SUPERLU_FREE(perm_c);
400  SUPERLU_FREE(R);
401  SUPERLU_FREE(C);
402  SUPERLU_FREE(etree);
403  if (lwork!=0) SUPERLU_FREE(work);
404  Destroy_CompCol_Matrix(&slu_A);
405  StatFree(&stat);
406  Destroy_SuperMatrix_Store(&slu_B);
407  Destroy_SuperMatrix_Store(&slu_X);
408  Destroy_SuperNode_Matrix(&slu_L);
409  Destroy_CompCol_Matrix(&slu_U);
410 #endif
411  }
412 
413  // check if DSAUPD failed
414  if (info<0)
415  {
416  if ((info<=-1)&&(info>=-6))
417  SG_SWARNING("ARPACK: DSAUPD failed. Wrong parameter passed.\n");
418  else if (info==-7)
419  SG_SWARNING("ARPACK: DSAUPD failed. Workaround array size is not sufficient.\n");
420  else
421  SG_SWARNING("ARPACK: DSAUPD failed. Error code: %d.\n", info);
422 
423  status = info;
424  }
425  else
426  {
427  if (info==1)
428  SG_SWARNING("ARPACK: Maximum number of iterations reached.\n");
429 
430  // allocate select for dseupd
431  int* select = SG_MALLOC(int, ncv);
432  // allocate d to hold eigenvalues
433  double* d = SG_MALLOC(double, 2*ncv);
434  // sigma for dseupd
435  double sigma = shift;
436 
437  // init ierr indicating dseupd possible errors
438  int ierr = 0;
439 
440  // specify that eigenvectors are going to be computed too
441  int rvec = 1;
442 
443  SG_SDEBUG("APRACK: Starting DSEUPD.\n");
444 
445  // call dseupd_ routine
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);
449 
450  // check for errors
451  if (ierr!=0)
452  {
453  SG_SWARNING("ARPACK: DSEUPD failed with status %d.\n", ierr);
454  status = -1;
455  }
456  else
457  {
458  SG_SDEBUG("ARPACK: Storing eigenpairs.\n");
459 
460  // store eigenpairs to specified arrays
461  for (i=0; i<nev; i++)
462  {
463  eigenvalues[i] = d[i];
464 
465  for (j=0; j<n; j++)
466  eigenvectors[j*nev+i] = v[i*n+j];
467  }
468  }
469 
470  // cleanup
471  SG_FREE(select);
472  SG_FREE(d);
473  }
474 
475  // cleanup
476  SG_FREE(all_);
477  SG_FREE(which_);
478  SG_FREE(resid);
479  SG_FREE(v);
480  SG_FREE(iparam);
481  SG_FREE(ipntr);
482  SG_FREE(workd);
483  SG_FREE(workl);
484 };
485 }
486 #endif /* HAVE_LAPACK */
487 #endif /* HAVE_ARPACK */

SHOGUN Machine Learning Toolbox - Documentation