SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros 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,
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 = 1;
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  StatInit(&stat);
207 
208  // factorize
209  slu_info = 0;
210  slu_Bv = doubleMalloc(n);
211  slu_Xv = doubleMalloc(n);
212  dCreate_Dense_Matrix(&slu_B,n,1,slu_Bv,n,SLU_DN,SLU_D,SLU_GE);
213  dCreate_Dense_Matrix(&slu_X,n,1,slu_Xv,n,SLU_DN,SLU_D,SLU_GE);
214  slu_B.ncol = 0;
215  SG_SDEBUG("SUPERLU: Factorizing\n");
216  dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
217  &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond, &ferr, &berr,
218  &mem_usage,&stat,&slu_info);
219  slu_B.ncol = 1;
220  if (slu_info)
221  {
222  SG_SERROR("SUPERLU: Failed to factorize matrix, got %d code\n", slu_info);
223  }
224  options.Fact = FACTORED;
225 #endif
226  }
227  else
228  {
229  // compute factorization according to pos value
230  if (pos)
231  {
232  // with Cholesky
233  SG_SDEBUG("ARPACK: Using Cholesky factorization.\n");
234  clapack_dpotrf(CblasColMajor,CblasUpper,n,matrix,n);
235  }
236  else
237  {
238  // with LUP
239  SG_SDEBUG("ARPACK: Using LUP factorization.\n");
240  ipiv = SG_CALLOC(int, n);
241  clapack_dgetrf(CblasColMajor,n,n,matrix,n,ipiv);
242  }
243  }
244  }
245  // main computation loop
246  SG_SDEBUG("ARPACK: Starting main computation DSAUPD loop.\n");
247  do
248  {
249  dsaupd_(&ido, bmat, &n, which_, &nev, &tol, resid,
250  &ncv, v, &ldv, iparam, ipntr, workd, workl,
251  &lworkl, &info);
252 
253  if ((ido==1)||(ido==-1)||(ido==2))
254  {
255  if (mode==1)
256  {
257  // compute (workd+ipntr[1]-1) = A*(workd+ipntr[0]-1)
258  cblas_dsymv(CblasColMajor,CblasUpper,
259  n,1.0,matrix,n,
260  (workd+ipntr[0]-1),1,
261  0.0,(workd+ipntr[1]-1),1);
262  }
263  if (mode==3)
264  {
265  if (!rhs)
266  {
267  if (use_superlu)
268  {
269 #ifdef HAVE_SUPERLU
270  // treat workd+ipntr(0) as B
271  for (i=0; i<n; i++)
272  slu_Bv[i] = (workd+ipntr[0]-1)[i];
273  slu_info = 0;
274  // solve
275  dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
276  &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond,
277  &ferr, &berr, &mem_usage, &stat, &slu_info);
278  if (slu_info)
279  SG_SERROR("SUPERLU: GOT %d\n", slu_info);
280  // move elements from resulting X to workd+ipntr(1)
281  for (i=0; i<n; i++)
282  (workd+ipntr[1]-1)[i] = slu_Xv[i];
283 #endif
284  }
285  else
286  {
287  for (i=0; i<n; i++)
288  (workd+ipntr[1]-1)[i] = (workd+ipntr[0]-1)[i];
289  if (pos)
290  clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
291  else
292  clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
293  }
294  }
295  else
296  // HAVE RHS
297  {
298  if (ido==-1)
299  {
300  if (use_superlu)
301  {
302 #ifdef HAVE_SUPERLU
303  for (i=0; i<n; i++)
304  slu_Bv[i] = (workd+ipntr[0]-1)[i];
305  slu_info = 0;
306  dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
307  &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond,
308  &ferr, &berr, &mem_usage, &stat, &slu_info);
309  if (slu_info)
310  SG_SERROR("SUPERLU: GOT %d\n", slu_info);
311  for (i=0; i<n; i++)
312  (workd+ipntr[1]-1)[i] = slu_Xv[i];
313 #endif
314  }
315  else
316  {
317  for (i=0; i<n; i++)
318  (workd+ipntr[1]-1)[i] = rhs[i]*(workd+ipntr[0]-1)[i];
319  if (pos)
320  clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
321  else
322  clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
323  }
324  }
325  if (ido==1)
326  {
327  if (use_superlu)
328  {
329 #ifdef HAVE_SUPERLU
330  for (i=0; i<n; i++)
331  slu_Bv[i] = (workd+ipntr[2]-1)[i];
332  slu_info = 0;
333  dgssvx(&options, &slu_A, perm_c, perm_r, etree, equed, R, C,
334  &slu_L, &slu_U, work, lwork, &slu_B, &slu_X, &rpg, &rcond,
335  &ferr, &berr, &mem_usage, &stat, &slu_info);
336  if (slu_info)
337  SG_SERROR("SUPERLU: GOT %d\n", slu_info);
338  for (i=0; i<n; i++)
339  (workd+ipntr[1]-1)[i] = slu_Xv[i];
340 #endif
341  }
342  else
343  {
344  for (i=0; i<n; i++)
345  (workd+ipntr[1]-1)[i] = (workd+ipntr[2]-1)[i];
346  if (pos)
347  clapack_dpotrs(CblasColMajor,CblasUpper,n,1,matrix,n,(workd+ipntr[1]-1),n);
348  else
349  clapack_dgetrs(CblasColMajor,CblasNoTrans,n,1,matrix,n,ipiv,(workd+ipntr[1]-1),n);
350  }
351  }
352  if (ido==2)
353  {
354  if (is_rhs_diag)
355  {
356  for (i=0; i<n; i++)
357  (workd+ipntr[1]-1)[i] = rhs[i]*(workd+ipntr[0]-1)[i];
358  }
359  else
360  {
361  cblas_dsymv(CblasColMajor,CblasUpper,
362  n,1.0,rhs,n,
363  (workd+ipntr[0]-1),1,
364  0.0,(workd+ipntr[1]-1),1);
365  }
366  }
367  }
368  }
369  }
370  } while ((ido==1)||(ido==-1)||(ido==2));
371 
372  if (!pos && mode==3) SG_FREE(ipiv);
373 
374  if (mode==3 && use_superlu)
375  {
376 #ifdef HAVE_SUPERLU
377  SUPERLU_FREE(slu_Bv);
378  SUPERLU_FREE(slu_Xv);
379  SUPERLU_FREE(perm_r);
380  SUPERLU_FREE(perm_c);
381  SUPERLU_FREE(R);
382  SUPERLU_FREE(C);
383  SUPERLU_FREE(etree);
384  if (lwork!=0) SUPERLU_FREE(work);
385  Destroy_CompCol_Matrix(&slu_A);
386  StatFree(&stat);
387  Destroy_SuperMatrix_Store(&slu_B);
388  Destroy_SuperMatrix_Store(&slu_X);
389  Destroy_SuperNode_Matrix(&slu_L);
390  Destroy_CompCol_Matrix(&slu_U);
391 #endif
392  }
393 
394  // check if DSAUPD failed
395  if (info<0)
396  {
397  if ((info<=-1)&&(info>=-6))
398  SG_SWARNING("ARPACK: DSAUPD failed. Wrong parameter passed.\n");
399  else if (info==-7)
400  SG_SWARNING("ARPACK: DSAUPD failed. Workaround array size is not sufficient.\n");
401  else
402  SG_SWARNING("ARPACK: DSAUPD failed. Error code: %d.\n", info);
403 
404  status = info;
405  }
406  else
407  {
408  if (info==1)
409  SG_SWARNING("ARPACK: Maximum number of iterations reached.\n");
410 
411  // allocate select for dseupd
412  int* select = SG_MALLOC(int, ncv);
413  // allocate d to hold eigenvalues
414  double* d = SG_MALLOC(double, 2*ncv);
415  // sigma for dseupd
416  double sigma = shift;
417 
418  // init ierr indicating dseupd possible errors
419  int ierr = 0;
420 
421  // specify that eigenvectors are going to be computed too
422  int rvec = 1;
423 
424  SG_SDEBUG("APRACK: Starting DSEUPD.\n");
425 
426  // call dseupd_ routine
427  dseupd_(&rvec, all_, select, d, v, &ldv, &sigma, bmat,
428  &n, which_, &nev, &tol, resid, &ncv, v, &ldv,
429  iparam, ipntr, workd, workl, &lworkl, &ierr);
430 
431  // check for errors
432  if (ierr!=0)
433  {
434  SG_SWARNING("ARPACK: DSEUPD failed with status %d.\n", ierr);
435  status = -1;
436  }
437  else
438  {
439  SG_SDEBUG("ARPACK: Storing eigenpairs.\n");
440 
441  // store eigenpairs to specified arrays
442  for (i=0; i<nev; i++)
443  {
444  eigenvalues[i] = d[i];
445 
446  for (j=0; j<n; j++)
447  eigenvectors[j*nev+i] = v[i*n+j];
448  }
449  }
450 
451  // cleanup
452  SG_FREE(select);
453  SG_FREE(d);
454  }
455 
456  // cleanup
457  SG_FREE(all_);
458  SG_FREE(which_);
459  SG_FREE(resid);
460  SG_FREE(v);
461  SG_FREE(iparam);
462  SG_FREE(ipntr);
463  SG_FREE(workd);
464  SG_FREE(workl);
465 };
466 }
467 #endif /* HAVE_LAPACK */
468 #endif /* HAVE_ARPACK */

SHOGUN Machine Learning Toolbox - Documentation