SCOREC / core

parallel finite element unstructured meshes
Other
183 stars 62 forks source link

apf::eigen not converging #131

Open zhang-alvin opened 6 years ago

zhang-alvin commented 6 years ago

I'm trying to get the eigenvectors and eigenvalues for a simple matrix using the apf::eigen() function:

  apf::Matrix3x3 testMat;
  testMat[0][0] = 1.0;
  testMat[0][1] = 0.0;
  testMat[0][2] = 15.0;
  testMat[1][0] = 0.0;
  testMat[1][1] = 1.0;
  testMat[1][2] = 0.0;
  testMat[2][0] = 0.0;
  testMat[2][1] = 2.0;
  testMat[2][2] = 5.0;

  apf::Vector3 eigenVectors[3];
  double eigenValues[3];
  apf::eigen(testMat, eigenVectors, eigenValues);

For this simple case, I would run into an assertion error at https://github.com/SCOREC/core/blob/master/apf/apfMatrix.cc#L76 because the algorithm had not converged. This seems abnormal since Octave finds the eigenvalues to be integers $\lambda = {1,5,1}$.

I compiled this code with the core libraries loaded via module on the SCOREC machines.

I also noticed that the maximum number of iterations is hard-coded as maxiters = 100 which seems like something that the user should decide.

zhang-alvin commented 6 years ago

I had forgotten about this issue. I remember talking to @mortezah about this and one limitation of apf::eigen is that the input matrix must be symmetric. The input might also have to be invertible. So I think what would be great is really a better error message that details such limitations.

AjinkyaDahale commented 5 years ago

Bringing this issue back from the dead a second time. In my specific case, which does have a symmetric matrix, changing maxiters = 100 to maxiters = 1000 does not change the matrix which is supposed to converge.

The culprit appears to be the convergence criterion used: https://github.com/SCOREC/core/blob/2ca5177fd172dd310a9412ad9b549b8afe291b2c/mth/mthQR.cc#L207-L218 In lines 211-212, the limit is set to an absolute value 1e-10. Instead the value should probably be relative to some value dependent on the matrix in question, here a. I'm looking for suggestions as to what that should be.

This does not work for the specific case @zhang-alvin mentioned, but that could just be because that matrix is not symmetric.