## Solvers not working

2 posts
• Page

**1**of**1**### Solvers not working

Hi,

The solver is meant to be part of an iterative algorithm where the (square) A and b in the Ax = b equations are slightly changed (only values, structure is exactly the same) between the iterations. The first iteration is solved without any problems, but during the second the solver will not provide an even remotely good answer (and no errors are reported).

I've reduced it to a minimal example, as the following (I'm also using the CUSP https://code.google.com/p/cusp-library/ library here to simplify handling of arrays and matrices, but I think you should be able to understand what's going on):

OS: latest Redhat

CUDA: 5.0

GPU: Telsa K20c

And here are the two A:s and b:s I'm using:

Iteration 1: (solves fine)

https://www.dropbox.com/s/84qf138lq8p240w/AAETT.mtx

https://www.dropbox.com/s/gbjvutqx21q3gw2/baffETT.mtx

Iteration 2: (does not solve)

https://www.dropbox.com/s/53638a6wtacl4h5/AATVA.mtx

https://www.dropbox.com/s/p1ixl70eo0m5xat/baffTVA.mtx

Does anyone have any idea what's going on?

The solver is meant to be part of an iterative algorithm where the (square) A and b in the Ax = b equations are slightly changed (only values, structure is exactly the same) between the iterations. The first iteration is solved without any problems, but during the second the solver will not provide an even remotely good answer (and no errors are reported).

I've reduced it to a minimal example, as the following (I'm also using the CUSP https://code.google.com/p/cusp-library/ library here to simplify handling of arrays and matrices, but I think you should be able to understand what's going on):

- Code: Select all
`#include "cula.hpp"`

#include <cusp/multiply.h>

#include <cusp/transpose.h>

#include <cusp/io/matrix_market.h>

#include <cusp/print.h>

using namespace std;

typedef double valueType;

typedef cusp::array1d<valueType, cusp::host_memory> Array1H;

typedef cusp::array1d<int, cusp::host_memory> iArray1H;

typedef cusp::array2d<valueType, cusp::host_memory> Array2H;

int main() {

culaInitialize();

Array2H A2_h;

Array1H A_h;

Array1H originalb_h;

cusp::io::read_matrix_market_file(A2_h, "AATVA.mtx");

const int n = A2_h.num_cols;

A_h = A2_h.values;

Array1H b_h(n);

Array1H temp_h(n);

cusp::io::read_matrix_market_file(originalb_h, "baffTVA.mtx");

cusp::blas::copy(originalb_h, b_h);

iArray1H ipiv_h(n);

culaStatus s = culaGesv(n, 1, &*A_h.begin(), n, &*ipiv_h.begin(), &*b_h.begin(), n); // b = A\b

if(s != culaNoError) {

if(s == culaArgumentError)

printf("Argument %d has an illegal value\n", culaGetErrorInfo());

else

printf("%s\n", culaGetStatusString(s));

}

cusp::multiply(A2_h, b_h, temp_h); // temp = A*b

cusp::blas::axpy(originalb_h, temp_h, -1.0f); // temp = temp - originalb

cusp::print(temp_h); // prints the error, that is b - A*x

culaShutdown();

return 0;

}

OS: latest Redhat

CUDA: 5.0

GPU: Telsa K20c

And here are the two A:s and b:s I'm using:

Iteration 1: (solves fine)

https://www.dropbox.com/s/84qf138lq8p240w/AAETT.mtx

https://www.dropbox.com/s/gbjvutqx21q3gw2/baffETT.mtx

Iteration 2: (does not solve)

https://www.dropbox.com/s/53638a6wtacl4h5/AATVA.mtx

https://www.dropbox.com/s/p1ixl70eo0m5xat/baffTVA.mtx

Does anyone have any idea what's going on?

- deveight
**Posts:**1**Joined:**Fri Jun 28, 2013 9:32 am

### Re: Solvers not working

These matrices have definitely been made dense, correct? And the symmetric portions have been populated? GESV assumes dense and fully filled, and the sparse MTX files are neither.

Keep in mind that GESV modifies the matrix and leaves the LU factorized matrix in its place, so if you call it twice with the same matrix it will not solve the second one properly. You need to use GETRS for subsequent solves against the same matrix.

Keep in mind that GESV modifies the matrix and leaves the LU factorized matrix in its place, so if you call it twice with the same matrix it will not solve the second one properly. You need to use GETRS for subsequent solves against the same matrix.

- john
- Administrator
**Posts:**587**Joined:**Thu Jul 23, 2009 2:31 pm

2 posts
• Page

**1**of**1**### Who is online

Users browsing this forum: No registered users and 1 guest