[Bf-blender-cvs] SVN commit: /data/svn/bf-blender [38650] branches/soc-2011-avocado/blender/ intern/autoseam: Successful calculation of eigen vectors using arpack.
shuvro sarker
shuvro05 at gmail.com
Sun Jul 24 07:24:46 CEST 2011
Revision: 38650
http://projects.blender.org/scm/viewvc.php?view=rev&root=bf-blender&revision=38650
Author: shuvro
Date: 2011-07-24 05:24:46 +0000 (Sun, 24 Jul 2011)
Log Message:
-----------
Successful calculation of eigen vectors using arpack. Now it is integrated with the present autoseam module. The only thing left is to alloacte the memory needed by the arpack functions dynamically. After doing that, the integration of arpack with autoseam will be complete.
Modified Paths:
--------------
branches/soc-2011-avocado/blender/intern/autoseam/AutoseamAdjacency.cpp
branches/soc-2011-avocado/blender/intern/autoseam/AutoseamAdjacency.h
Modified: branches/soc-2011-avocado/blender/intern/autoseam/AutoseamAdjacency.cpp
===================================================================
--- branches/soc-2011-avocado/blender/intern/autoseam/AutoseamAdjacency.cpp 2011-07-24 04:34:46 UTC (rev 38649)
+++ branches/soc-2011-avocado/blender/intern/autoseam/AutoseamAdjacency.cpp 2011-07-24 05:24:46 UTC (rev 38650)
@@ -31,6 +31,8 @@
#include "BKE_report.h"
#include "AutoseamAdjacency.h"
+//#define EIGEN_DEBUG_ARPACK 1
+
#ifdef WITH_ARPACK
#include "dssimp.h"
//#include "f2c.h"
@@ -39,7 +41,7 @@
#define THRESHOLD_ZERO 0.0005
-#define NUM_EIGEN_VAL 30
+//#define NUM_EIGEN_VAL 25
#define minimum(a,b) (a <= b) ? a : b;
@@ -168,7 +170,7 @@
* The function will calculate eigen values using the functions
* of arpack.
*/
-int AutoseamAdjacency:: calculate_eigen_arpack()
+double * AutoseamAdjacency:: calculate_eigen_arpack()
{
#ifdef WITH_ARPACK
@@ -218,22 +220,26 @@
static cilist io___35 = { 0, 6, 0, 0, 0 };
- nx = matrix_dimension;
- //nx = 10;
- n = nx * nx; // here n is the matrix size
+ // This is the dimension of the matrix
+ n = matrix_dimension;
-
- // The number of eigen values we want to calculate
- // This should be the minimum number between matrix_dimension and NUM_EIGEN_VAL
- //nev = 4;
- nev = minimum(matrix_dimension, NUM_EIGEN_VAL);
- ncv = 20;
+ /**
+ *The number of eigen values we want to calculate
+ * This should be the minimum number between matrix_dimension and NUM_EIGEN_VAL
+ * The constraints - nev <= maxnev
+ * nev + 1 <= ncv <= maxncv
+ * n >= ncv
+ */
+ nev = minimum(matrix_dimension - 1, NUM_EIGEN_VAL);
+ num_eigen_values = nev;
+ ncv = minimum(nev + 1, 20);
+
*(unsigned char *)bmat = 'I';
- //we are looking for the largest eigen values
- s_copy(which, "LM", (ftnlen)2, (ftnlen)2);
+ //we are looking for the smallest eigen values
+ s_copy(which, "SM", (ftnlen)2, (ftnlen)2);
if (n > 256) {
s_wsle(&io___7);
@@ -241,7 +247,7 @@
e_wsle();
goto L9000;
}
- else if (nev > 10) {
+ else if (nev > NUM_EIGEN_VAL) {
s_wsle(&io___8);
do_lio(&c__9, &c__1, " ERROR with _SSIMP: NEV is greater than MAXNEV ", (ftnlen)47);
e_wsle();
@@ -283,7 +289,7 @@
if (ido == -1 || ido == 1) {
//perform matrix-vector multiplication
//av_(&nx, &workd[ipntr[0] - 1], &workd[ipntr[1] - 1]);
- mul_matrix_vector(nx, &workd[ipntr[0] - 1], &workd[ipntr[1] - 1]);
+ mul_matrix_vector(n, &workd[ipntr[0] - 1], &workd[ipntr[1] - 1]);
goto L10;
}
@@ -309,45 +315,25 @@
}
else {
- // no major error
- rvec = TRUE_;
+ // no major error
+ rvec = TRUE_;
- dseupd_(&rvec, "All", select, d__, v, &c__256, &sigma, bmat, &n,
+ dseupd_(&rvec, "All", select, d__, v, &c__256, &sigma, bmat, &n,
which, &nev, &tol, resid, &ncv, v, &c__256, iparam, ipntr,
workd, workl, &lworkl, &ierr, (ftnlen)3, (ftnlen)1, (ftnlen)2);
-
- // Block to print eigen values and eigen vectors
-
- #ifdef EIGEN_DEBUG_ARPACK
- {
- int i,j;
-
- for(i = 1; i <= matrix_dimension; i++)
- {
- printf("eigen value : %f\n", d__[i]);
- //now the eigen vector
- printf("The eigen vector %d :\n", i);
-
- for(j = 0; j < matrix_dimension; j++){
- printf("%lf ", *(v + i + (25*j) ) );
- }
- }
- }
- #endif
-
- // Now we are going to print the eigen values and eigen vectors here.
- /* %----------------------------------------------% */
- /* | Eigenvalues are returned in the first column | */
- /* | of the two dimensional array D and the | */
- /* | corresponding eigenvectors are returned in | */
- /* | the first NCONV (=IPARAM(5)) columns of the | */
- /* | two dimensional array V if requested. | */
- /* | Otherwise, an orthogonal basis for the | */
- /* | invariant subspace corresponding to the | */
- /* | eigenvalues in D is returned in V. | */
- /* %----------------------------------------------% */
+ /* %----------------------------------------------% */
+ /* | Eigenvalues are returned in the first column | */
+ /* | of the two dimensional array D and the | */
+ /* | corresponding eigenvectors are returned in | */
+ /* | the first NCONV (=IPARAM(5)) columns of the | */
+ /* | two dimensional array V if requested. | */
+ /* | Otherwise, an orthogonal basis for the | */
+ /* | invariant subspace corresponding to the | */
+ /* | eigenvalues in D is returned in V. | */
+ /* %----------------------------------------------% */
+
if (ierr != 0) {
s_wsle(&io___32);
@@ -366,23 +352,32 @@
}
else {
-
- nconv = iparam[4];
- i__1 = nconv;
- for (j = 1; j <= i__1; ++j) {
+
+ // Block to print eigen values and eigen vectors
+
+ #ifdef EIGEN_DEBUG_ARPACK
+ {
+ int i,j;
- //matrix-vector multiplication
- //av_(&nx, &v[(j << 8) - 256], ax);
- mul_matrix_vector(nx, &v[(j << 8) - 256], ax);
- d__1 = -d__[j - 1];
-
- daxpy_(&n, &d__1, &v[(j << 8) - 256], &c__1, ax, &c__1);
- d__[j + 24] = dnrm2_(&n, ax, &c__1);
- d__[j + 24] /= (d__1 = d__[j - 1], abs(d__1));
-
- /* L20: */
+ for(i = 0; i < nev; i++)
+ {
+ printf("eigen value : %f\n", d__[i]);
+ printf("The eigen vector %d :\n", i);
+
+ //prints one eigen vector
+ for(j = 0; j < matrix_dimension; j++){
+ printf("%lf ", v[i*256 + j]);
+ }
+ }
}
+ #endif
+ //set the eigen values
+ int i;
+ eigen_values = Eigen::VectorXd(nev);
+ for(i = 0; i < nev; i++) eigen_values(i) = d__[i];
+ // return the eigen vectors
+ return v;
}
@@ -394,7 +389,7 @@
#endif
L9000:
- return 0;
+ return NULL;
@@ -403,12 +398,14 @@
bool AutoseamAdjacency::generate_seam()
{
+ double *eigen_vectors;
MatrixXd& a = m_adjacency;
clear_eigenspaces();
int n = a.rows();
+
for (int i=0; i<n;++i) {
a(i,i) = 0.0;
double rowsum = a.row(i).sum();
@@ -417,26 +414,39 @@
// We need to calculate eigen values from here.
- calculate_eigen_arpack();
+ eigen_vectors = calculate_eigen_arpack();
//dssimp_();
- if(a.rows() && a.cols()){
+ if(a.rows() && a.cols() && eigen_vectors != NULL){
- Eigen::SelfAdjointEigenSolver<MatrixXd> es(a);
- Eigen::VectorXd evalues(es.eigenvalues());
+ //Eigen::SelfAdjointEigenSolver<MatrixXd> es(a);
+ //Eigen::VectorXd evalues(es.eigenvalues());
int f = 0;
bool found = false;
MatrixXd aplus;
MatrixXd aminus;
- while ( (f < n)) {
+ while ( (f < num_eigen_values)) {
// Eigenvalues seem to be sorted largest to smallest, we need the 30 smallest
// in the future only those will be calculated by the algorithm (if we use ARPACK)
- if(fabs(evalues[n-f-1]) > 0.0005){
+
+ //if(fabs(eigen_values[n-f-1]) > 0.0005){
+ if(fabs(eigen_values[f]) > 0.0005){
+
+ int loop;
+ //now construct a vector with the values
+ Eigen::VectorXd eigen_vector = Eigen::VectorXd(matrix_dimension);
+
+ for(loop = 0; loop < matrix_dimension; loop++){
+ // This hardcoded value will be removed later
+ eigen_vector(loop) = *(eigen_vectors + f*256 + loop);
+ }
+
- AutoseamEigenspace* aes = new AutoseamEigenspace(evalues[n-f-1], es.eigenvectors().col(n-f-1));
+ //AutoseamEigenspace* aes = new AutoseamEigenspace(evalues[n-f-1], es.eigenvectors().col(n-f-1));
+ AutoseamEigenspace* aes = new AutoseamEigenspace(eigen_values[f], eigen_vector);
// split Eigenspace into two subspaces F+ and F- where the ith entry in the eigenvector is positive/negative
aes->split();
aes->fill_adjacency(a, aplus, aminus);
Modified: branches/soc-2011-avocado/blender/intern/autoseam/AutoseamAdjacency.h
===================================================================
--- branches/soc-2011-avocado/blender/intern/autoseam/AutoseamAdjacency.h 2011-07-24 04:34:46 UTC (rev 38649)
+++ branches/soc-2011-avocado/blender/intern/autoseam/AutoseamAdjacency.h 2011-07-24 05:24:46 UTC (rev 38650)
@@ -33,7 +33,8 @@
#include "AutoseamEigenspace.h"
#include <vector>
-#define NUM_EIGEN_VECTOR 30
+#define NUM_EIGEN_VAL 25
+//#define NUM_EIGEN_VECTOR 30
class AutoseamAdjacency
{
@@ -57,7 +58,7 @@
float get_value(int row, int col);
void debug(std::ofstream& fout);
- int calculate_eigen_arpack();
+ double * calculate_eigen_arpack();
int mul_matrix_vector(int, double *, double *);
@@ -71,8 +72,12 @@
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list