Page 1 of 1

### Solvers not working

Posted: Wed Jul 03, 2013 6:23 am
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):
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;

const int n = A2_h.num_cols;
A_h = A2_h.values;

Array1H b_h(n);
Array1H temp_h(n);
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?

### Re: Solvers not working

Posted: Wed Jul 03, 2013 10:02 am
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.