anderkve / FYS3150

https://anderkve.github.io/FYS3150
26 stars 14 forks source link

spsolve() is our enemy #109

Closed siljevik closed 1 year ago

siljevik commented 1 year ago

Hi, When we are trying to use spsolve(), we get this message:

libc++abi: terminating with uncaught exception of type std::runtime_error: spsolve(): solution not found zsh: abort ./main.exe

And we have tried to google a lot, and we don't understand Armadillos pages either.. It does not show how stuff work, it should be more examples.

We use it this way: arma::cx_vec u_np1 = spsolve(spA, b);

where 'spA' is a 'sp_cx_mat' and 'b' is a 'cx_vec'.

anderkve commented 1 year ago

Hi!

Are you sure the problem is an issue with armadillo, and not a problem with how spA and b are constructed? I.e., could it be that there's a problem with spA or b in such a way that the equation Ax=b simply has no solution, and the message from armadillo is correct?

As an example, the following test code gives the same error, because the equation I try to solve simply has no solution:

#include <armadillo>
#include <iostream>

using namespace std;
using namespace arma;

int main()
{
  sp_cx_mat A = sp_cx_mat(2, 2);
  A(0,0) = 1.0;
  A(1,1) = 0.0;

  cx_vec b = {1.0, 2.0};

  cx_vec x = spsolve(A,b);

  cout << "x = " << endl;
  cout << x << endl;

  return 0;
}

On my Ubuntu machine running this code gives the (correct) error

terminate called after throwing an instance of 'std::runtime_error'
  what():  spsolve(): solution not found
Aborted (core dumped)
siljevik commented 1 year ago

If I print spA, it looks like this: (0, 0) (nan,nan) (1, 0) (nan,nan) (18, 0) (nan,nan) (0, 1) (nan,nan) (1, 1) (nan,nan) (2, 1) (nan,nan) (19, 1) (nan,nan) (1, 2) (nan,nan) (2, 2) (nan,nan) (20, 2) (nan,nan) (3, 3) (nan,nan) (4, 3) (nan,nan) (21, 3) (nan,nan) (3, 4) (nan,nan) (4, 4) (nan,nan) (5, 4) (nan,nan) (22, 4) (nan,nan) (4, 5) (nan,nan) (5, 5) (nan,nan) (23, 5) (nan,nan) (6, 6) (nan,nan) (7, 6) (nan,nan) (24, 6) (nan,nan) etc..

So the A is what is wrong here then? How do we insert values to an sp_cx_mat? We did it like for an cx_mat, but maybe it is wrong? i.e. spA(k,k) = some value or A.diag(M-1) = some other value

We made our code first using normal matrices, then we changed to cx_mat, and then we figured out we needed to use the sp_cx_mat. Maybe a little late here, and that may be what's causing all this trouble.

anderkve commented 1 year ago

Yeah, if your matrix is full of nan ("not a number") values, then using it in a matrix equation (or anything else) is bound to fail.

It should be perfectly fine to fill sp_cx_mat objects with the syntax you have used.

I would rather check if the values you fill A with are nan for some reason. For instance, if you have a problem in the computation of the parameter "r" (see project description) then that could affect every element of the A matrix you construct.

siljevik commented 1 year ago

Yes, that is where our issue is! But, I don't understand how to use i in the equation?? r is an cx_double, like this: arma::cx_double r = (sqrt(-1))(dt)/(2.0(h*h));

Is there somehow another way we should use i (sqrt(-1))?

anderkve commented 1 year ago

Try this:

    cx_double i = cx_double(0,1);
    cx_double r = (i * dt) / (2.0 * h * h);

or simply this

   cx_double r = (1.0i * dt) / (2.0 * h * h);
anderkve commented 1 year ago

For the latter option, you may have to compile with the flag -std=c++14.

siljevik commented 1 year ago

Yay! It works now!:D Thank you!