PLearn 0.1
parpack.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // PLearn (A C++ Machine Learning Library)
00004 // Copyright (C) 1998 Pascal Vincent
00005 // Copyright (C) 1999-2002 Pascal Vincent, Yoshua Bengio and University of Montreal
00006 //
00007 
00008 // Redistribution and use in source and binary forms, with or without
00009 // modification, are permitted provided that the following conditions are met:
00010 // 
00011 //  1. Redistributions of source code must retain the above copyright
00012 //     notice, this list of conditions and the following disclaimer.
00013 // 
00014 //  2. Redistributions in binary form must reproduce the above copyright
00015 //     notice, this list of conditions and the following disclaimer in the
00016 //     documentation and/or other materials provided with the distribution.
00017 // 
00018 //  3. The name of the authors may not be used to endorse or promote
00019 //     products derived from this software without specific prior written
00020 //     permission.
00021 // 
00022 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00023 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00024 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00025 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00027 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00028 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00029 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00030 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00031 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00032 // 
00033 // This file is part of the PLearn library. For more information on the PLearn
00034 // library, go to the PLearn Web site at www.plearn.org
00035 
00036 
00037 #ifndef parpack_INC
00038 #define parpack_INC
00039 
00040 #include "arpack_proto.h"
00041 #include "Mat.h"
00042 
00043 namespace PLearn {
00044 using namespace std;
00045 
00046 
00088 template<class MatT>
00089 int eigenSparseSymmMat(MatT& A, Vec& e_values, Mat& e_vectors, FORTRAN_Integer& n_evalues,
00090                        int max_n_iter=300, bool compute_vectors=true, bool largest_evalues=true,
00091                        bool according_to_magnitude=true, bool both_ends=false, real ncv2nev_ratio=1.5)
00092 {
00093 #ifdef NOARPACK
00094     PLERROR("eigenSparseSymmMat: ARPACK not available on this system!");
00095     return 0;
00096 #else
00097     FORTRAN_Integer ido=0;
00098     char bmat[1];
00099     bmat[0] = 'I';
00100     char which[2];
00101     FORTRAN_Integer n=A.length();
00102     if (e_vectors.length()!=n_evalues || e_vectors.width()!=n)
00103         PLERROR("eigenSparseSymmMat: expected e_vectors.width(%d)=A.length(%d), e_vectors.length(%d)=e_values.length(%d)",
00104                 e_vectors.width(),n,e_vectors.length(),n_evalues);
00105     if (both_ends)
00106         strncpy(which,"BE",2); 
00107     else
00108     {
00109         which[0]= largest_evalues? 'L' : 'S';
00110         which[1]= according_to_magnitude? 'M' : 'A';
00111     }
00112     real tol=0;
00113     FORTRAN_Integer ncv=MIN(3+int(n_evalues*ncv2nev_ratio),n-1);
00114     e_vectors.resize(ncv,n); 
00115     FORTRAN_Integer iparam[11];
00116     iparam[0]=1;
00117     iparam[2]=max_n_iter;
00118     iparam[6]=1;
00119     FORTRAN_Integer ipntr[11];
00120     Vec workd(3*n);
00121     FORTRAN_Integer lworkl = (ncv * ncv) + (ncv * 8);
00122     Vec workl(lworkl);
00123     Vec resid(n);
00124     FORTRAN_Integer info=0;
00125     for (;;) {
00126 #ifdef USEDOUBLE
00127         dsaupd_(&ido, bmat, &n, which, &n_evalues, &tol, resid.data(), &ncv, e_vectors.data(), &n,
00128                 iparam, ipntr, workd.data(), workl.data(), &lworkl, &info, 1, 2);
00129 #endif
00130 #ifdef USEFLOAT
00131         ssaupd_(&ido, bmat, &n, which, &n_evalues, &tol, resid.data(), &ncv, e_vectors.data(), &n,
00132                 iparam, ipntr, workd.data(), workl.data(), &lworkl, &info, 1, 2);
00133 #endif
00134         if (ido == -1 || ido == 1) {
00135             Vec x=workd.subVec(ipntr[0]-1,n);
00136             Vec z=workd.subVec(ipntr[1]-1,n);
00137             product(A, x, z);
00138         } else break;
00139     }
00140     if (info != 0 && info != 1)
00141     {
00142         PLWARNING("eigenSparseSymmMat: saupd returning error %ld",info);
00143         return info;
00144     }
00145     if (info > 0)
00146     {
00147         n_evalues = iparam[4];
00148         e_values.resize(n_evalues);
00149     }
00150     e_vectors.resize(n_evalues,n);
00151     if (n_evalues>0)
00152     {
00153         FORTRAN_Integer rvec = compute_vectors;
00154         TVec<FORTRAN_Integer> select(ncv);
00155         FORTRAN_Integer ierr;
00156         real sigma =0;
00157 #ifdef USEDOUBLE
00158         dseupd_(&rvec, "All", select.data(), e_values.data(), e_vectors.data(), &n, &sigma, 
00159                 bmat, &n, which, &n_evalues, &tol, resid.data(), &ncv, e_vectors.data(), &n, 
00160                 iparam, ipntr, workd.data(), workl.data(), &lworkl, &ierr, 3, 1, 2);
00161 #endif
00162 #ifdef USEFLOAT
00163         sseupd_(&rvec, "All", select.data(), e_values.data(), e_vectors.data(), &n, &sigma, 
00164                 bmat, &n, which, &n_evalues, &tol, resid.data(), &ncv, e_vectors.data(), &n, 
00165                 iparam, ipntr, workd.data(), workl.data(), &lworkl, &ierr, 3, 1, 2);
00166 #endif
00167         if (ierr != 0)
00168         {
00169             PLWARNING("eigenSparseSymmMat: seupd returning error %ld",ierr);
00170             return ierr;
00171         }
00172     }
00173 #endif
00174     return info;
00175 }
00176 
00177 
00185 template<class MatT>
00186 int eigenSparseNonSymmMat(MatT& A, Vec e_values, Mat e_vectors, FORTRAN_Integer& n_evalues,
00187                           int max_n_iter=300, bool compute_vectors=true, bool largest_evalues=true,
00188                           bool according_to_magnitude=true, bool both_ends=false)
00189 {
00190 #ifdef NOARPACK
00191     PLERROR("eigenSparseNonSymmMat: ARPACK not available on this system!");
00192     return 0;
00193 #else
00194     FORTRAN_Integer ido=0;
00195     char bmat[1];
00196     bmat[0] = 'I';
00197     char which[2];
00198     FORTRAN_Integer n=A.length();
00199     if (e_vectors.length()!=n_evalues || e_vectors.width()!=n)
00200         PLERROR("eigenSparseNonSymmMat: expected e_vectors.width(%d)=A.length(%d), e_vectors.length(%d)=e_values.length(%d)",
00201                 e_vectors.width(),n,e_vectors.length(),n_evalues);
00202     if (both_ends)
00203         strncpy(which,"BE",2); 
00204     else
00205     {
00206         which[0]= largest_evalues? 'L' : 'S';
00207         which[1]= according_to_magnitude? 'M' : 'R';//according to magnitude or according to real part
00208     }
00209     real tol=0;
00210     FORTRAN_Integer ncv=MIN(3+int(n_evalues*1.5),n-1);
00211     e_vectors.resize(ncv,n); 
00212     FORTRAN_Integer iparam[11];
00213     iparam[0]=1;
00214     iparam[2]=max_n_iter;
00215     iparam[6]=1;
00216     FORTRAN_Integer ipntr[11];
00217     Vec workd(3*n);
00218     FORTRAN_Integer lworkl = 3*(ncv * ncv) + (ncv * 6);
00219     Vec workl(lworkl);
00220     Vec resid(n);
00221     FORTRAN_Integer info=0;
00222     for (;;) {
00223 #ifdef USEDOUBLE
00224         dnaupd_(&ido, bmat, &n, which, &n_evalues, &tol, resid.data(), &ncv, e_vectors.data(), &n,
00225                 iparam, ipntr, workd.data(), workl.data(), &lworkl, &info, 1, 2);
00226 #endif
00227 #ifdef USEFLOAT
00228         snaupd_(&ido, bmat, &n, which, &n_evalues, &tol, resid.data(), &ncv, e_vectors.data(), &n,
00229                 iparam, ipntr, workd.data(), workl.data(), &lworkl, &info, 1, 2);
00230 #endif
00231         if (ido == -1 || ido == 1) {
00232             Vec x=workd.subVec(ipntr[0]-1,n);
00233             Vec z=workd.subVec(ipntr[1]-1,n);
00234             product(A, x, z);
00235         } else break;
00236     }
00237     if (info != 0 && info != 1)
00238     {
00239         PLWARNING("eigenSparseNonSymmMat: naupd returning error %ld",info);
00240         return info;
00241     }
00242     Vec e_valuesIm(n_evalues+1);
00243     Vec workev(3*ncv);
00244     if (info > 0)
00245     {
00246         n_evalues = iparam[4];
00247         e_values.resize(n_evalues+1);
00248     }
00249     e_vectors.resize(n_evalues+1,n);
00250     if (n_evalues>0)
00251     {
00252         FORTRAN_Integer rvec = compute_vectors;
00253         TVec<FORTRAN_Integer> select(ncv);
00254         FORTRAN_Integer ierr;
00255         real sigmai =0;
00256         real sigmar =0;
00257 #ifdef USEDOUBLE
00258         dneupd_(&rvec, "A", select.data(), e_values.data(), e_valuesIm.data(), e_vectors.data(), &n,
00259                 &sigmar, &sigmai, workev.data(), bmat, &n, which, &n_evalues, &tol, 
00260                 resid.data(), &ncv, e_vectors.data(), &n, iparam, ipntr, workd.data(), 
00261                 workl.data(), &lworkl, &ierr, 3, 1, 2);
00262 #endif
00263 #ifdef USEFLOAT
00264         sneupd_(&rvec, "A", select.data(), e_values.data(), e_valuesIm.data(), e_vectors.data(), &n,
00265                 &sigmar, &sigmai, workev.data(), bmat, &n, which, &n_evalues, &tol, 
00266                 resid.data(), &ncv, e_vectors.data(), &n, iparam, ipntr, workd.data(), 
00267                 workl.data(), &lworkl, &ierr, 3, 1, 2);
00268 #endif
00269         if (ierr != 0)
00270         {
00271             PLWARNING("eigenSparseNonSymmMat: neupd returning error %ld",ierr);
00272             return ierr;
00273         }
00274     }
00275 #endif
00276     return info;
00277 }
00278 
00279 
00280 } // end of namespace PLearn
00281 
00282 
00283 #endif
00284 
00285 
00286 /*
00287   Local Variables:
00288   mode:c++
00289   c-basic-offset:4
00290   c-file-style:"stroustrup"
00291   c-file-offsets:((innamespace . 0)(inline-open . 0))
00292   indent-tabs-mode:nil
00293   fill-column:79
00294   End:
00295 */
00296 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines