| sp_auxlib_meat.hpp | | sp_auxlib_meat.hpp | |
|---|
| | | | |
| skipping to change at line 52 | | skipping to change at line 52 | |
|---|
| podarray<blas_int> iparam, ipntr; | | podarray<blas_int> iparam, ipntr; | |
| podarray<eT> rwork; // Not used in this case. | | podarray<eT> rwork; // Not used in this case. | |
| | | | |
| run_aupd(n_eigvals, p, true /* sym, not gen */, n, tol, resid, ncv, v,
ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | | run_aupd(n_eigvals, p, true /* sym, not gen */, n, tol, resid, ncv, v,
ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | |
| | | | |
| if(info != 0) | | if(info != 0) | |
| { | | { | |
| return false; | | return false; | |
| } | | } | |
| | | | |
|
| // The process has converged, and now we need to recover the actual | | // The process has converged, and now we need to recover the actual eig | |
| // eigenvectors using seupd(). | | envectors using seupd(). | |
| blas_int rvec = 1; // .TRUE | | blas_int rvec = 1; // .TRUE | |
| blas_int nev = n_eigvals; | | blas_int nev = n_eigvals; | |
| char howmny = 'A'; | | char howmny = 'A'; | |
| char bmat = 'I'; // We are considering the standard eigenvalue problem. | | char bmat = 'I'; // We are considering the standard eigenvalue problem. | |
| char which[3] = "LM"; // We want the eigenvalues with the largest magni
tude. | | char which[3] = "LM"; // We want the eigenvalues with the largest magni
tude. | |
| podarray<blas_int> select(ncv); // Logical array of dimension NCV. | | podarray<blas_int> select(ncv); // Logical array of dimension NCV. | |
| blas_int ldz = n; | | blas_int ldz = n; | |
| | | | |
| // seupd() will output directly into the eigval and eigvec objects. | | // seupd() will output directly into the eigval and eigvec objects. | |
| eigval.set_size(n_eigvals); | | eigval.set_size(n_eigvals); | |
| eigvec.set_size(n, n_eigvals); | | eigvec.set_size(n, n_eigvals); | |
| | | | |
| arpack::seupd(&rvec, &howmny, select.memptr(), eigval.memptr(), eigvec.
memptr(), &ldz, (eT*) NULL, &bmat, &n, which, &nev, &tol, resid.memptr(), &
ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), wor
kl.memptr(), &lworkl, &info); | | arpack::seupd(&rvec, &howmny, select.memptr(), eigval.memptr(), eigvec.
memptr(), &ldz, (eT*) NULL, &bmat, &n, which, &nev, &tol, resid.memptr(), &
ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd.memptr(), wor
kl.memptr(), &lworkl, &info); | |
| | | | |
| // Check for errors. | | // Check for errors. | |
| if(info != 0) | | if(info != 0) | |
| { | | { | |
|
| | | std::stringstream tmp; | |
| | | tmp << "eigs_sym(): ARPACK error " << info << " in seupd()"; | |
| | | arma_debug_warn(true, tmp.str()); | |
| return false; | | return false; | |
|
| // std::stringstream tmp; | | | |
| // tmp << "eigs_sym(): ARPACK error " << info << " in seupd()"; | | | |
| // arma_stop(tmp.str()); | | | |
| } | | } | |
| | | | |
| return (info == 0); | | return (info == 0); | |
| } | | } | |
| #else | | #else | |
| { | | { | |
| arma_ignore(eigval); | | arma_ignore(eigval); | |
| arma_ignore(eigvec); | | arma_ignore(eigvec); | |
| arma_ignore(X); | | arma_ignore(X); | |
| arma_ignore(n_eigvals); | | arma_ignore(n_eigvals); | |
| | | | |
| skipping to change at line 132 | | skipping to change at line 131 | |
|---|
| podarray<blas_int> iparam, ipntr; | | podarray<blas_int> iparam, ipntr; | |
| podarray<T> rwork; // Not used in the real case. | | podarray<T> rwork; // Not used in the real case. | |
| | | | |
| run_aupd(n_eigvals, p, false /* gen, not sym */, n, tol, resid, ncv, v,
ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | | run_aupd(n_eigvals, p, false /* gen, not sym */, n, tol, resid, ncv, v,
ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | |
| | | | |
| if(info != 0) | | if(info != 0) | |
| { | | { | |
| return false; | | return false; | |
| } | | } | |
| | | | |
|
| // The process has converged, and now we need to recover the actual | | // The process has converged, and now we need to recover the actual eig | |
| // eigenvectors using neupd(). | | envectors using neupd(). | |
| blas_int rvec = 1; // .TRUE | | blas_int rvec = 1; // .TRUE | |
| blas_int nev = n_eigvals; | | blas_int nev = n_eigvals; | |
| char howmny = 'A'; | | char howmny = 'A'; | |
| char bmat = 'I'; // We are considering the standard eigenvalue problem. | | char bmat = 'I'; // We are considering the standard eigenvalue problem. | |
| char which[3] = "LM"; // We want the eigenvalues with the largest magni
tude. | | char which[3] = "LM"; // We want the eigenvalues with the largest magni
tude. | |
| podarray<blas_int> select(ncv); // Logical array of dimension NCV. | | podarray<blas_int> select(ncv); // Logical array of dimension NCV. | |
| podarray<T> dr(nev + 1); // Real array of dimension NEV + 1. | | podarray<T> dr(nev + 1); // Real array of dimension NEV + 1. | |
| podarray<T> di(nev + 1); // Real array of dimension NEV + 1. | | podarray<T> di(nev + 1); // Real array of dimension NEV + 1. | |
| podarray<T> z(n * (nev + 1)); // Real N by NEV array if HOWMNY = 'A'. | | podarray<T> z(n * (nev + 1)); // Real N by NEV array if HOWMNY = 'A'. | |
| blas_int ldz = n; | | blas_int ldz = n; | |
| podarray<T> workev(3 * ncv); | | podarray<T> workev(3 * ncv); | |
| | | | |
| arpack::neupd(&rvec, &howmny, select.memptr(), dr.memptr(), di.memptr()
, z.memptr(), &ldz, (T*) NULL, (T*) NULL, workev.memptr(), &bmat, &n, which
, &nev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipnt
r.memptr(), workd.memptr(), workl.memptr(), &lworkl, rwork.memptr(), &info)
; | | arpack::neupd(&rvec, &howmny, select.memptr(), dr.memptr(), di.memptr()
, z.memptr(), &ldz, (T*) NULL, (T*) NULL, workev.memptr(), &bmat, &n, which
, &nev, &tol, resid.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipnt
r.memptr(), workd.memptr(), workl.memptr(), &lworkl, rwork.memptr(), &info)
; | |
| | | | |
| // Check for errors. | | // Check for errors. | |
| if(info != 0) | | if(info != 0) | |
| { | | { | |
|
| | | std::stringstream tmp; | |
| | | tmp << "eigs_gen(): ARPACK error " << info << " in neupd()"; | |
| | | arma_debug_warn(true, tmp.str()); | |
| return false; | | return false; | |
|
| // std::stringstream tmp; | | | |
| // tmp << "eigs_gen(): ARPACK error " << info << " in neupd()"; | | | |
| // arma_stop(tmp.str()); | | | |
| } | | } | |
| | | | |
| // Put it into the outputs. | | // Put it into the outputs. | |
| eigval.set_size(n_eigvals); | | eigval.set_size(n_eigvals); | |
| eigvec.set_size(n, n_eigvals); | | eigvec.set_size(n, n_eigvals); | |
| | | | |
| for (uword i = 0; i < n_eigvals; ++i) | | for (uword i = 0; i < n_eigvals; ++i) | |
| { | | { | |
| eigval[i] = std::complex<T>(dr[i], di[i]); | | eigval[i] = std::complex<T>(dr[i], di[i]); | |
| } | | } | |
| | | | |
| skipping to change at line 255 | | skipping to change at line 253 | |
|---|
| podarray<blas_int> iparam, ipntr; | | podarray<blas_int> iparam, ipntr; | |
| podarray<T> rwork; | | podarray<T> rwork; | |
| | | | |
| run_aupd(n_eigvals, p, false /* gen, not sym */, n, tol, resid, ncv, v,
ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | | run_aupd(n_eigvals, p, false /* gen, not sym */, n, tol, resid, ncv, v,
ldv, iparam, ipntr, workd, workl, lworkl, rwork, info); | |
| | | | |
| if(info != 0) | | if(info != 0) | |
| { | | { | |
| return false; | | return false; | |
| } | | } | |
| | | | |
|
| // The process has converged, and now we need to recover the actual | | // The process has converged, and now we need to recover the actual eig | |
| // eigenvectors using neupd(). | | envectors using neupd(). | |
| blas_int rvec = 1; // .TRUE | | blas_int rvec = 1; // .TRUE | |
| blas_int nev = n_eigvals; | | blas_int nev = n_eigvals; | |
| char howmny = 'A'; | | char howmny = 'A'; | |
| char bmat = 'I'; // We are considering the standard eigenvalue problem. | | char bmat = 'I'; // We are considering the standard eigenvalue problem. | |
| char which[3] = "LM"; // We want the eigenvalues with the largest magni
tude. | | char which[3] = "LM"; // We want the eigenvalues with the largest magni
tude. | |
| podarray<blas_int> select(ncv); // Logical array of dimension NCV. | | podarray<blas_int> select(ncv); // Logical array of dimension NCV. | |
| podarray<std::complex<T> > d(nev + 1); // Real array of dimension NEV +
1. | | podarray<std::complex<T> > d(nev + 1); // Real array of dimension NEV +
1. | |
| podarray<std::complex<T> > z(n * nev); // Real N by NEV array if HOWMNY
= 'A'. | | podarray<std::complex<T> > z(n * nev); // Real N by NEV array if HOWMNY
= 'A'. | |
| blas_int ldz = n; | | blas_int ldz = n; | |
| podarray<std::complex<T> > workev(2 * ncv); | | podarray<std::complex<T> > workev(2 * ncv); | |
| | | | |
| skipping to change at line 279 | | skipping to change at line 276 | |
|---|
| eigval.set_size(n_eigvals); | | eigval.set_size(n_eigvals); | |
| eigvec.set_size(n, n_eigvals); | | eigvec.set_size(n, n_eigvals); | |
| std::complex<T> sigma; | | std::complex<T> sigma; | |
| | | | |
| arpack::neupd(&rvec, &howmny, select.memptr(), eigval.memptr(), | | arpack::neupd(&rvec, &howmny, select.memptr(), eigval.memptr(), | |
| (std::complex<T>*) NULL, eigvec.memptr(), &ldz, (std::complex<T>*) &sigma,
(std::complex<T>*) NULL, workev.memptr(), &bmat, &n, which, &nev, &tol, res
id.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd
.memptr(), workl.memptr(), &lworkl, rwork.memptr(), &info); | | (std::complex<T>*) NULL, eigvec.memptr(), &ldz, (std::complex<T>*) &sigma,
(std::complex<T>*) NULL, workev.memptr(), &bmat, &n, which, &nev, &tol, res
id.memptr(), &ncv, v.memptr(), &ldv, iparam.memptr(), ipntr.memptr(), workd
.memptr(), workl.memptr(), &lworkl, rwork.memptr(), &info); | |
| | | | |
| // Check for errors. | | // Check for errors. | |
| if(info != 0) | | if(info != 0) | |
| { | | { | |
|
| | | std::stringstream tmp; | |
| | | tmp << "eigs_gen(): ARPACK error " << info << " in neupd()"; | |
| | | arma_debug_warn(true, tmp.str()); | |
| return false; | | return false; | |
|
| // std::stringstream tmp; | | | |
| // tmp << "eigs_gen(): ARPACK error " << info << " in neupd()"; | | | |
| // arma_stop(tmp.str()); | | | |
| } | | } | |
| | | | |
| return (info == 0); | | return (info == 0); | |
| } | | } | |
| #else | | #else | |
| { | | { | |
| arma_ignore(eigval); | | arma_ignore(eigval); | |
| arma_ignore(eigvec); | | arma_ignore(eigvec); | |
| arma_ignore(X); | | arma_ignore(X); | |
| arma_ignore(n_eigvals); | | arma_ignore(n_eigvals); | |
| | | | |
| skipping to change at line 340 | | skipping to change at line 337 | |
|---|
| // "NCV must satisfy the two inequalities 2 <= NCV-NEV and NCV <= N". | | // "NCV must satisfy the two inequalities 2 <= NCV-NEV and NCV <= N". | |
| // "It is recommended that NCV >= 2 * NEV". | | // "It is recommended that NCV >= 2 * NEV". | |
| ncv = (2 * nev < n) ? 2 * nev : ((nev + 2 < n) ? nev + 2 : n); | | ncv = (2 * nev < n) ? 2 * nev : ((nev + 2 < n) ? nev + 2 : n); | |
| v.set_size(n * ncv); // Array N by NCV (output). | | v.set_size(n * ncv); // Array N by NCV (output). | |
| rwork.set_size(ncv); // Work array of size NCV for complex calls. | | rwork.set_size(ncv); // Work array of size NCV for complex calls. | |
| ldv = n; // "Leading dimension of V exactly as declared in the calling
program." | | ldv = n; // "Leading dimension of V exactly as declared in the calling
program." | |
| | | | |
| // IPARAM: integer array of length 11. | | // IPARAM: integer array of length 11. | |
| iparam.zeros(11); | | iparam.zeros(11); | |
| iparam(0) = 1; // Exact shifts (not provided by us). | | iparam(0) = 1; // Exact shifts (not provided by us). | |
|
| iparam(2) = 1000; // Maximum iterations; all the examples use 300, but | | iparam(2) = 1000; // Maximum iterations; all the examples use 300, but | |
| they | | they were written in the ancient times. | |
| // were written in the ancient times. | | | |
| iparam(6) = 1; // Mode 1: A * x = lambda * x. | | iparam(6) = 1; // Mode 1: A * x = lambda * x. | |
| | | | |
| // IPNTR: integer array of length 14 (output). | | // IPNTR: integer array of length 14 (output). | |
| ipntr.set_size(14); | | ipntr.set_size(14); | |
| | | | |
| // Real work array used in the basic Arnoldi iteration for reverse | | // Real work array used in the basic Arnoldi iteration for reverse | |
| // communication. | | // communication. | |
| workd.set_size(3 * n); | | workd.set_size(3 * n); | |
| | | | |
| // lworkl must be at least 3 * NCV^2 + 6 * NCV. | | // lworkl must be at least 3 * NCV^2 + 6 * NCV. | |
| | | | |
| skipping to change at line 412 | | skipping to change at line 408 | |
|---|
| // Nothing to do here, things have converged. | | // Nothing to do here, things have converged. | |
| break; | | break; | |
| default: | | default: | |
| { | | { | |
| return; // Parent frame can look at the value of info. | | return; // Parent frame can look at the value of info. | |
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
| // The process has ended; check the return code. | | // The process has ended; check the return code. | |
|
| if (info != 0 && info != 1) | | if( (info != 0) && (info != 1) ) | |
| { | | { | |
|
| | | // Print warnings if there was a failure. | |
| | | std::stringstream tmp; | |
| | | | |
| | | if(sym) | |
| | | { | |
| | | tmp << "eigs_sym(): ARPACK error " << info << " in saupd()"; | |
| | | } | |
| | | else | |
| | | { | |
| | | tmp << "eigs_gen(): ARPACK error " << info << " in naupd()"; | |
| | | } | |
| | | | |
| | | arma_debug_warn(true, tmp.str()); | |
| | | | |
| return; // Parent frame can look at the value of info. | | return; // Parent frame can look at the value of info. | |
| } | | } | |
| } | | } | |
| #endif | | #endif | |
| } | | } | |
| | | | |
| End of changes. 12 change blocks. |
|---|
| 19 lines changed or deleted | | 32 lines changed or added | |
|---|
|