PLearn 0.1
|
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 :