Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

How to solve Ax=b , where A is a symmetry matrix . #28

Open
ztdepztdep opened this issue Dec 2, 2021 · 9 comments
Open

How to solve Ax=b , where A is a symmetry matrix . #28

ztdepztdep opened this issue Dec 2, 2021 · 9 comments

Comments

@ztdepztdep
Copy link

Which function should i use.

@michael-lehn
Copy link
Owner

If A is positiv definite use the cholesky factorization (flens::lapack::potrf). Lapack documentation: http://www.netlib.org/lapack/explore-html/d1/d7a/group__double_p_ocomputational_ga2f55f604a6003d03b5cd4a0adcfb74d6.html

If A is not positive definite use the Bunch-Kaufman diagonal pivoting method (flens::sy::trf). Lapack documentation:
https://www.netlib.org/lapack/explore-html/d3/db6/group__double_s_ycomputational_gad91bde1212277b3e909eb6af7f64858a.html

@ztdepztdep
Copy link
Author

If A is positiv definite use the cholesky factorization (flens::lapack::potrf). Lapack documentation: http://www.netlib.org/lapack/explore-html/d1/d7a/group__double_p_ocomputational_ga2f55f604a6003d03b5cd4a0adcfb74d6.html

If A is not positive definite use the Bunch-Kaufman diagonal pivoting method (flens::sy::trf). Lapack documentation: https://www.netlib.org/lapack/explore-html/d3/db6/group__double_s_ycomputational_gad91bde1212277b3e909eb6af7f64858a.html

Thank you very much. In my problem, I need to solve Ax=b repeatly where the matrix A doesn't change.
So how to seperate the solve process into two stages. first decomposet A, then solve it with different right hand side.

@michael-lehn
Copy link
Owner

It is already separated. flens::lapack::potrf and flens::sy::trf are computing the factorization. flens::lapack::potrs and flens::sy::trs solve the equation and require that the matrix was previously factorized.

Exactly how this is handled by LAPACK with 1) dgepotrf and dgepotrs or 2) dgesytrf and dgesytrs:

  1. http://www.netlib.org/lapack/explore-html/d1/d7a/group__double_p_ocomputational_ga2f55f604a6003d03b5cd4a0adcfb74d6.html

and

http://www.netlib.org/lapack/explore-html/d1/d7a/group__double_p_ocomputational_ga167aa0166c4ce726385f65e4ab05e7c1.html

  1. https://www.netlib.org/lapack/explore-html/d3/db6/group__double_s_ycomputational_gad91bde1212277b3e909eb6af7f64858a.html

and

https://www.netlib.org/lapack/explore-html/d3/db6/group__double_s_ycomputational_ga6a223e61effac7232e98b422f147f032.html

@ztdepztdep
Copy link
Author

Thanks.
If I solve Ax=b with the "lapack::sv(A, piv, b)" , is the solution time slower than the " flens::sy::trf" and "flens::sy::trs". I am working on a spectral multidomain code, the speed is import for me to improve the efficiency.

@michael-lehn
Copy link
Owner

If I solve Ax=b with the "lapack::sv(A, piv, b)" , is the solution time slower than the " flens::sy::trf" and "flens::sy::trs".

Is this a question or an observation? What type does A have? GeMatrix or SyMatrix?

@ztdepztdep
Copy link
Author

ztdepztdep commented Dec 5, 2021

I tired to compile the "lapack-potrs.cc", the compiler feedback the following error:
`./flens/storage/array/arrayview.tcc:200:30: error: invalid use of member function ‘flens::ArrayView<T, I, A>::IndexType flens::ArrayView<T, I, A>::firstIndex() const [with T = double; I = flens::IndexOptions<>; A = std::allocator; flens::ArrayView<T, I, A>::IndexType = int]’ (did you forget the ‘()’ ?)
_data += _firstIndex - firstIndex;

/flens/storage/array/arrayview.tcc:202:17: error: cannot convert ‘flens::ArrayView<double, flens::IndexOptions<>, std::allocator >::firstIndex’ from type ‘flens::ArrayView<double, flens::IndexOptions<>, std::allocator >::IndexType (flens::ArrayView<double, flens::IndexOptions<>, std::allocator >::)() const {aka int (flens::ArrayView<double, flens::IndexOptions<>, std::allocator >::)() const}’ to type ‘flens::ArrayView<double, flens::IndexOptions<>, std::allocator >::IndexType {aka int}’
_firstIndex = firstIndex;`

Could you please help me out.

@michael-lehn
Copy link
Owner

There seemed to be an issue with gcc 11. After fixing an missing include I can compile lapack-potrs.cc in flens/examples as follows:

g++-11 -Wall -std=c++11 -I /Users/lehn/work/FLENS lapack-potrs.cc

This is not using an external LAPACK.

  1. Can you first confirm that this compiles on you machine? If not: what compiler are you using
  2. For using an external LAPACK: What LAPACK implementation are you using? MKL?

@ztdepztdep
Copy link
Author

I compiled it in the KDevelop 4 with cmake, it gives me the above error. I still cann't make it work.

I tried to compile it in the terminal wth :
g++ -std=c++11 -O3 -I /home/ztdep/Downloads/FLENS-2012-10-01 main.cpp -o b it feeds back correct results.

@michael-lehn
Copy link
Owner

I am not familiar with KDevelop 4. But I compiling it on the terminal works it is just a matter of configuring the IDE. Most likely you just have to specify the correct compiler flags somewhere.

Is there a way to see how KDevelop 4 invokes the compiler?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants