SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
arpack.cpp
Go to the documentation of this file.
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 */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation